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1.  OVERVIEW 


INTRODUCTION 

As  the  title  suggests,  the  work  described  in  this  report  is  motivated  by  the  problem 
of  geolocating  the  transmitter  or  receiver  of  a  satellite-relayed  signal.  The  situation  is 
that  pictured  in  Figure  1.1;  the  doppler  shift  between  a  transmitter  and  receiver  is  meas¬ 
ured  over  a  period  of  time  and  utilized  to  determine  the  receiver  or  transmitter  location. 
A  number  of  other  parameters,  in  addition  to  the  unknown  location  of  the  transmitter, 
are  intimately  involved  in  the  computation.  Among  these  are  the  transmitter  frequency, 
the  satellite’s  orbital  parameters,  course  and  speed  of  moving  stations,  etc.  Such  param¬ 
eters  are  not  fundamentally  different  from  position  variables,  and,  in  fact,  any  of  the 
entire  set  of  parameters  may  be  considered  part  of  the  unknown  state  to  be  estimated.1 
The  remaining  parameters  must  be  known  a  priori,  and  their  accuracy  will  affect  that  of 
the  state  estimate. 

Although  many  of  the  details  of  the  present  study  are  specific  to  the  above  prob¬ 
lem,  our  methodology  and  results  are  quite  general.  Consequently,  we  begin  with  a  brief 
description  of  the  abstract  estimation  problem,  illustrating  it  by  specializing  to  geoloca¬ 
tion.  Let  x  be  a  vector  whose  components  xs,  s  =  1,...,  N  represent  a  set  of  parameters 
to  be  estimated.  We  suppose  the  existence  of  additional  parameters  q  which  together 
with  the  state  vector  x  determine  (in  the  absence  of  noise)  the  measurements.  This  rela¬ 
tionship  is  described  by  a  set  of  functions,  mj(x,  q),  i  —  1,  ...,  M,  where  M  is  the  number 
of  measurements.  The  actual  measurements  m;  will  be  noisy,  resulting  in  errors 
Anij  =  mj(x,  q)  -  m-v  We  choose  as  an  estimate  a  state  xt  which  minimizes  these  errors. 
More  precisely,  xe  is  the  nonlinear  least  squares  estimate  determined  by  minimizing  the 
cost  function 


cm  -  sfofr-  . 

i  <7i 


(1-1) 


C(xe)  =  minC(x'' . 


(1.2) 


The  weights  erf  are  the  variances  of  the  measurement  noise.2  In  the  above  geolocation 
problem,  the  state  x  would  be  the  transmitter  latitude  and  longitude,  and  m;  would  be 
the  doppler-shifted  frequency  at  the  receiver  at  time  tj,  which  we  write  as  f(x,  q,  tj).  The 
weights  crj2  correspond  to  the  variances  of  the  measured  receiver  frequencies  1';.  It  is 
sometimes  also  desirable  to  include  other  measurement  types,  such  as  satellite  azimuth 
and  elevation,  which  may  be  used  to  simultaneously  determine  the  orbit  and  geolocate. 

Let  us  consider  measures  of  the  quality  of  the  estimate  obtained  by  (1.2).  The  com¬ 
ponents  of  xt  are  random  variables  since  they  are  functions  of  the  random  variables  m ;. 
Thus,  an  appropriate  error  measure  for  a  component  (xe)s  of  xe  is  its  mean  squared  error, 

1  The  measurements  may  not  always  support  a  solution,  however. 

2  This  particular  formulation  is  appropriate  for  measurements  which  are  stochastically  independent, 
a  restriction  which  is  relaxed  in  Section  2.  If  the  noise  is  also  Gaussian,  then  (l.l)-(1.2)  is  a  max¬ 
imum  likelihood  estimator. 
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i.e,  E  (xs-(xe)s  )2,  where  E  denotes  expectation  and  x  is  the  true  state.  More  generally, 
we  may  utilize  the  covariance  matrix 


xAE((*-x,)(*-x«)t) 


(1.3) 


as  a  measure  of  error.  (The  dagger  f  denotes  the  transpose  of  the  column  vector  x.) 

While  (1.3)  reflects  the  errors  due  to  measurement  noise,  it  is  also  important  to 
determine  the  sensitivity  of  xe  to  errors  in  the  parameters  q.  This  relationship  will  gen¬ 
erally  be  deterministic  (but,  see  (2.28)  -  (2.36)).  A  given  error  ^  in  the  kth  parameter 
will  produce  an  error  6x  in  the  state  estimate.  To  the  first  order,  this  dependence  is 
linear  so  that  the  total  error  maj.  be  expressed  in  the  form 


<5x  =  G  <5q  , 


(1.4) 


where  G  is  a  matrix  determined  by  the  true  values,  i.e.,  G  =  G(x,  q).  Typical  com¬ 
ponents  of  q  for  the  geolocation  problem  might  include  satellite  orbit  inclination, 
transmitter  frequency,  etc. 

The  parallel  between  the  problem  statement  and  the  sensitivity/error  analysis  is 
illustrated  in  Figure  1.2.  A  given  problem  is  realized  by  specifying  which  variables  are  to 
be  estimated,  and  which  are  inputs.  For  the  sensitivity  we  also  need  estimates  of  the 
accuracies  of  the  inputs.  Note  the  parenthetical  quantities  in  the  figure.  For  the  estima¬ 
tion  problem,  the  measurement  variances  may  be  used  to  weight  the  influence  of  the 
inputs  on  the  estimate  (cf.  equation  (l.l)).  For  the  sensitivity,  the  values  of  the  parame¬ 
ters  and  measurements  are  used  to  determine  the  relationship  of  the  estimate  to  inputs3 
and,  thereby,  the  transfer  function  between  the  input  accuracies  and  the  estimate  accu¬ 
racies. 

In  the  remainder  of  the  report,  a  weak  attempt  has  been  made  to  place  the  sections 
in  order  of  increasing  mathematical  detail.  The  current  section,  although  specific  to  the 
geolocation  problem,  contains  a  minimum  of  mathematical  analysis.  The  subsection  fol¬ 
lowing  this  introduction  specifies  the  measurement  functions,  i.e.,  the  dependency  of  the 
doppler  on  the  various  parameters  x  and  q.  The  third  subsection  contains  a  qualitative 
discussion  of  sensitivity  in  the  same  context. 

Section  2  is  a  self-contained  development  of  the  sensitivity  and  error  analysis  for 
the  general  estimation  problem.  Part  of  that  section  also  describes  an  effective  numerical 
technique  for  the  minimization  (1.2).  Solving  these  problems  requires  the  computation  of 
the  measurement  functions  and  their  derivatives.  Such  computations  are  detailed  in 
Section  3  for  a  stationary  transmitter  and  receiver  on  an  oblate  earth.  This  model  is 
extended  in  Section  4  to  include  ground  station  motion.  Section  4  (accompanied  by 
Appendix  A)  may  be  read  independently  of  the  rest  of  this  report.  Computations 

3  The  values  are  needed  (as  well  as  their  accuracies)  because  the  sensitivity  is  generally  a  nonlinear 
function  of  the  state  and/or  the  parameters  If  the  true  state  is  known,  as  in  a  simulation,  then  the 
measurement  values  are  not  needed. 


3 


(meas  variances) 

parameters- 
measurements- 


( parameters ,  state  or  measurements) 

parameter  accuracies- 
measurement  variances- 


estimate  accuracy 

Figure  1.2.  Generic  diagram  of  a  state  estimator  and  sensitivity  analysis.  Note  that  the 
sensitivity  for  an  unknown  state  is  approximated  by  using  the  state  estimate  obtained 
from  the  measurements. 


f  A  nonlinear  least  squares  estimator  in  the  current  context. 
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requiring  a  specific  satellite  orbit  model  are  relegated  to  Appendix  B,  where  they  are 
detailed  for  Keplerian  motion.  The  actual  software  can  accept  an  arbitrary  orbit  model, 
in  which  case  it  computes  the  derivatives  numerically. 


DOPPLER  MEASUREMENTS 

Let  us  return  to  Figure  1.1.  The  arrows  marked  ”  doppler  up”  and  ’’doppler  down” 
are  intended  to  indicate  the  doppler  shift  in  signal  frequency  due  to  the  relative  motion 
of  the  satellite  and  the  transmitter  or  receiver,  respectively.  For  example,  if  we  denote 
the  transmitter  frequency  by  f0  and  the  relative  doppler  shift  on  the  uplink  by  by,  the 
frequency  of  the  signal  received  by  the  satellite  will  be  (ignoring  relativistic  effects) 


f  =  f0  (1  +  bu)  . 


(1.5) 


The  relevant  variables  for  the  calculation  of  b  are  found  in  Figure  1.3.  The  vector 
R  represents  the  position  of  the  ground  station  relative  to  the  center  of  the  earth,  r  is 
that  of  the  satellite,  and 


pir-R 


(1.6) 


is  the  vector  between  them.  The  relative  doppler  between  a  source  and  a  receiver  is  pro¬ 
portional  to  the  time  derivative  of  the  distance  between  them.  More  precisely,  we  have 

b=-ii||p||,  (1.7) 


where  c  is  the  speed  of  light.  Denoting  the  time  derivative  and  scalar  product  by  dots, 
we  have,  from  (1.7), 


b  — i  4r(P-P),/!  =  ---r£r-P 

c  dt  c  |  p  I 


Note  that  equation  (1.8)  is  coordinate-system-independent.4  From  (1.6) 


p  =  r  -  R  . 


(1.8) 


(1.9) 


Equations  (1.8)  and  (1.9)  give  b  as  a  function  of  station  motion  (R,  R)  and  satellite 

4  In  contrast,  the  following  equation,  (1.9),  depends  or.  the  coordinate  system  in  the  sense  that  R 
and  r  do  (i.e.,  a  moving  coordinate  system  will  affect  the  velocities).  Of  course  their  difference,  p, 
does  not. 
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motion  (r,  r).  Typical  computations  take  place  in  the  geocentric  coordinate  system  (so- 
called  IJK  inertial  coordinates  [l]).  Consequently,  even  if  a  station  is  stationary  with 
respect  to  the  earth,  R  is  not  zero  because  of  the  earth’s  rotation.  Then,  given  the 
station’s  latitude,  longitude,  and  altitude  at  some  fixed  point  in  time,  as  well  as  addi¬ 
tional  parameters  such  as  course  and  speed  if  there  is  ground  motion,  we  can  compute 
R(t)  and  R(t)  at  an  arbitrary  time  t.  As  mentioned  in  the  previous  subsection,  any  of 
these  parameters  may  be  considered  either  as  part  of  the  unknown  state  x  or  as  given 
values  q  in  the  measurement  functions  m;(x,  q)  of  (1.1). 

As  with  R,  the  satellite  position  and  velocity,  r(t)  and  r(t),  may  be  computed  util¬ 
izing  a  finite  set  of  parameters  modeling  the  satellite’s  motion.  For  example,  the  Kepler 
model,  detailed  in  Appendix  B,  contains  six  parameters:  the  orbit  semimajor  axis,  its 
eccentricity,  the  mean  anomaly,  inclination,  longitude  of  the  ascending  node,  and  the 
argument  of  perigee.  The  epoch  (i.e.,  point  in  time)  at  which  these  hold  must  also  be 
known.  For  example,  in  the  geolocation  problem,  one  might  take  the  unknown  state  x 
to  be  the  2-dimensional  vector  whose  components  are  transmitter  latitude  and  longitude, 
and  q  to  comprise  12  parameters,  the  transmitter  altitude,  receiver  latitude,  longitude, 
and  altitude,  the  transmitter  frequency  f0,  the  6  Kepler  elements,  and  a  satellite  relay 
(offset)  frequency  which  we  denote  fso  (in  all  12  dimensions). 

The  above  discussion  has  outlined  the  functional  dependencies  bup(x,  q)  and 
bDown(xi  <0  f°r  the  uplink  and  downlink  relative  doppler  shifts.  It  is  now  a  short  step  to 
connect  these  to  fR,  the  frequency  determined  at  the  receiver.  The  signal  leaves  the 
transmitter  with  a  frequency  f0  and  arrives  at  the  satellite  with  a  frequency  f0  (1  +  by). 
This  signal  is  then  relayed  by  the  satellite  with  a  possible  offset  of  fso,  resulting  in  an 
emitted  signal  of  frequency  f0  (l  +  by)  -f-  fSO  This  propagates  and,  due  to  the  downlink 
doppler,  arrives  at  the  receiver  with  a  frequency  [f0  (1  +  by)  +  fs0][l  +  bD].  Finally,  the 
receiver  introduces  another  offset,  which  we  write  -fR0  (usually  set  close  in  value  to 
-  (fo  +  fso)i  >n  order  that  the  receiver  output  frequency  be  as  close  as  possible  to  the 
total  doppler  shift).  In  summary,  the  output  frequency  of  the  receiver  is 

=  [fo(f  +  M  +  fso]  I1  +  bD]  -  fR0  .  (l-!0) 

We  have  thus  arrived  at  a  point  at  which  we  can  construct  the  measurement  function(s) 
fR(x,  q,  tj)  by  substituting  the  dependencies  bu(x,  q,  tj)  and  bR(x,  q,  tj)  into  (1.10).  The 
residuals  m;(x,  q)  -  mj  of  (1.1)  are  given  by  fR(x,  q,  tj)  -  f  j,  where  the  f;  are  the  frequen¬ 
cies  measured  at  the  receiver  output  during  an  experiment. 

There  is  substantial  value  in  writing  (1.10)  in  a  slightly  different  form.  First,  we 
rearrange  the  terms 


fR  =  b(jfo(l  +  bp)  +  bo(fo  +  fso)  +  (fso  +  fo  -  fRo)  •  (f -f  f) 


As  mentioned  above,  fRo  is  chosen  to  make  the  last  term  close  to  zero.  Regardless,  all 
the  information  lies  in  the  size  of  the  doppler  shift,  which  is  of  the  order  of  b  f0.  An 
error  in  f0  causes  a  substantial  percentage  error  in  this  quantity.  In  fact,  since  b  (by  and 
bp  are  of  the  same  order)  is  much  less  than  1,  the  percentage  error  is  approximately 


1  % 


<rf /bf0  —  may  remedy  this  by  estimating  most  of  the  uncertainty  in  the 

b  f0 

base  frequency  f0  (as  well  as  that  of  fso  and  fR0)  through  an  additional  state  variable, 


fB  =  fso  +  fo  -  fR0  • 


(1.12) 


Equation  (1.11)  becomes 


1r  =  bu  f0(l  4-  bp)  +  bp  (f0  4-  fso)  +  1b  •  (1-13) 


Since  fR  is  now  a  state  variable,  the  percentage  error  in  fR  of  equation  (1.13)  is  approxi- 
2<7f()b  <T(0  j  <7f() 

mately  — — —  «  which  is  much  less  than  - — — .  Of  course,  the  price  we  pay  is  the 

biQ  1q  b  Iq 

estimation  of  the  additional  state  variable  fB. 

Finally,  we  may  write  (1.13)  more  conveniently  by  dividing  by  f0  to  get 


d  Ail 

“  fo 

fpA 

==  bu  (1  +  bD)  +  bD  (1  +  -j— )  4-  ge  ,  (1*14) 

*o 


where  the  (possibly)  unknown  parameter  gB  is 


A  fB 

Kb  =  7~ 
•o 


fso  _  fRO 
fo  fo 


and  the  corresponding  measurements  are 


(1.15) 


(1.16) 


Even  if  f0  is  unknown,  but  a  rough  value  is  available,  equations  (1.14)  and  (1.16)  will 
often  be  a  very  accurate  approximation,  with  fo  as  the  rough  a  priori  value  and  gR 
estimating  the  true  value  of  f0  via  equation  (1.15). 

The  cost  function  corresponding  to  the  above  measurements  is 


C(x)  = 


J  \2 


(d(x,  q,  tj)  -  d ,) 

1  „  2 

^d, 


(1.17) 
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where,  from  (1.16),  the  variance  of  dj  is  cr^  =  The  minimum  of  the  cost  function 

(1.17)  satisfies  the  set  of  N  nonlinear  equations  VxC(xe)  =  0  and  may  be  solved  by  any 
number  of  steepest  descent  procedures.  We  recommend  the  technique  described  at  the 
end  of  section  2  (equations  (2.37)  to  (2.40)),  which  has  proved  highly  effective  for  the 
present  problem. 


SENSITIVITY  /  ERRORS 

In  the  geolocation  problem,  as  in  almost  all  estimation  problems,  it  is  critical  to 
have  an  an  idea  of  the  accuracy  of  the  solution  obtained.  An  unusually  poor  geometry 
may  produce  an  ill-conditioned  set  of  equations  whose  solution  is  of  essentially  no  practi¬ 
cal  value  (e.g.,  the  transmitter  is  in  Detroit  with  an  average  error  of  ±30°  latitude  and 
±80°  longitude).  Even  in  less  serious  cases,  extremely  erroneous  conclusions  may  be 
drawn  from  using  a  location  estimate  without  taking  into  account  its  accuracy.  Further¬ 
more,  the  geolocation  accuracy  is  not  a  single  fixed  value,  but  varies  with  the  satellite’s 
orbit,  transmitter  and  receiver  locations,  measurement  errors,  accuracy  of  the  orbit 
model,  etc.  These  explicit  dependencies  are  almost  always  relevant. 

Extensive  software  was  developed  to  determine  the  sensitivity  to  input  parameters 
and  provide  estimates  of  the  geolocation  accuracy.  The  methodology,  which  is  fully 
described  in  Section  2,  estimates  the  solution  errors  analytically  without  a  need  for 
extensive  Monte  Carlo  runs.  In  fact,  if  the  geolocation  estimates  are  derived  by  a 
steepest  descent  solution  of  the  nonlinear  least  squares  problem,  the  sensitivities  may  be 
provided  at  negligible  additional  computation.  Regardless,  independent  of  the  minimiza¬ 
tion  technique,  the  sensitivity  computation  is  relatively  fast  and  suitable  for  operation  in 
real  time. 

The  inputs  and  outputs  of  the  sensitivity  are  diagramed  in  Figures  1.4  through  1.6. 
The  desired  output  is  the  expected  error  (actually  the  RMS  error)  in  the  transmitter 
position  estimate.  In  Figure  1.4  this  is  depicted  as  a  function  of  the  orbit  errors,  the 
errors  in  frequency  parameters,  the  transmitter  location  error,  and  the  receiver  measure¬ 
ment  errors.  In  many  cases  it  is  also  necessary  to  determine  or  redetermine  the  satellite 
orbit  elements.  Such  a  determination  represents  an  inversion  of  the  original  problem: 
the  satellite  elements  (parameters)  become  the  unknown  state,  and  the  transmitter  loca¬ 
tion  is  included  as  a  known  input  parameter  qk.  In  that  context  (Figure  1.5),  the  sensi¬ 
tivity  software  outputs  the  accuracies  of  the  estimates  of  the  orbit  elements  which,  in 
turn,  may  be  used  as  inputs  to  the  sensitivity  analysis  of  the  geolocation  problem  (Figure 
1.4).  Of  course,  this  procedure  presupposes  the  existence  of  a  second  transmitter  whose 
location  is  known  and  may  thus  be  used  in  the  determination  of  the  orbit  elements.  It  is 
also  theoretically  possible  to  simultaneously  determine  the  transmitter  and  orbit  parame¬ 
ters,  thereby  avoiding  the  need  for  a  second  transmitter  (Figure  1.6);  however,  this 
requires  a  relatively  amicable  geometry  and/or  good  azimuth  and  elevation  data. 

We  now  describe  in  somewhat  more  detail  the  operation  of  the  sensitivity  software. 
Table  1.1  contains  a  listing  of  the  relevant  variables.  As  briefly  indicated  in  the  intro¬ 
duction,  we  distinguish  between  what  we  term  a  parameter  q,  which  has  a  fixed  error  for 

6  Satellite  azimuth  and  elevation  relative  to  a  given  ground  station  represent  additional  measure¬ 
ments  (mj  =  rn(x,  q,  tj))  which  are  useful  when  the  orbit  is  part  of  the  state  X, 
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Figure  1.4.  Inputs  and  ouputs  of  the  sensitivity  for  the  problem  of  determining  receiver  location 
and/or  motion. 
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Figure  1.5.  Inputs  and  ouputs  of  the  sensitivity  for  the  problem  of  determining  satellite  orbit 
elements. 
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Figure  1.6.  Inputs  and  ouputs  of  the  sensitivity  for  the  problem  of  determining  receiver  location 
and  satellite  orbit  elements. 
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Table  1.1.  Sensitivity  Parameters 


Variable 


Type 


transmitter 

latitude,  longitude 
altitude 
course,  speed 
frequency 


parameters 

parameter 

parameters 

parameter 


receiver 

latitude,  longitude,  altitude 
course,  speed 
frequency(t;)t 


parameters 

parameters 

measurements 


satellite 

elements' ' 

offset  frequency 
azimuth(tj),  elevation (tiU 
position(t|),  velocity(tj)  **  T 


parameters 
parameter 
measu  remen  ts 
(stochastic)  parameters 


f  Here,  t;  is  the  time  of  the  ith  frequency  measurement  fj.  Actually,  tj  is  also  a  measured 
quantity,  so  that  there  are  two  measurements,  f;  and  tj,  for  each  i. 

ft  mean  motion,  eccentricity,  mean  anomaly,  inclination,  ascending  node,  perigee,  etc. 

fff  Observe  that  these  represent  six  numbers  at  each  time  tj  since  position  and  velocity 
are  vectors. 
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any  one  geolocation  problem,  and  a  measurement  variable,  which  is  a  random  variable 
representing  a  set  of  noisy  measurements  m;,  whose  errors  share  a  common  probability 
distribution.  An  example  of  the  former  is  the  inclination  of  the  satellite  orbit,  while  the 
set  of  receiver  frequency  measurements  falls  in  the  latter  category.  More  precisely,  meas¬ 
urements  are  those  variables  whose  residuals  are  minimized  in  a  cost  function  in  order  to 
estimate  the  state  (cf.,  equation  (1.1)).  Sometimer  the  errors  of  a  set  of  parameters  will 
also  have  a  common  distribution;  we  term  these  ’’stochastic  parameters.”6  The  useful¬ 
ness  of  this  concept  is  discussed  below. 

Now  that  we  have  identified  the  variables  and  their  roles,  let  us  be  a  little  more 
specific  about  what  we  mean  by  their  ’’accuracies.”  The  error  in  a  parameter  is  con¬ 
sidered  to  be  fixed  throughout  an  experiment.  Thus,  if  the  transmitter  frequency  Is 
recorded  as  10  Mhz  when  it  is  actually  9.999  Mhz,  then  the  parameter  error  is  0.001 
Mhz.  Corresponding  to  such  an  error,  the  sensitivity  analysis  will  supply  an  error  for 
the  output  estimati ;  for  example,  a  latitude  error  of  0.02°.  If  there  are  no  measurement 
errors  (i.e.,  just  parameter  errors),  then  the  output  error  is  deterministic  and  a  linear 
function  (to  the  first  order)  of  the  input  error.  In  the  present  example,  an  error  of  -0.002 
Mhz  in  transmitter  frequency  would  produce  a  latitude  error  of  -0.04°. 

On  the  other  hand,  the  measurement  errors  are  stochastic.  Thus,  we  might  say 
that  the  set  of  received  frequencies  has  a  root  mean  square  (RMS)  error  of  erf  =  3Hz, 
meaning  that  if  we  average  the  squares  of  the  errors  of  the  measurements  and  take  the 
square  root  of  that  average,  we  get  3  Hz.  Correspondingly,  the  output  estimate  will 
have  a  random  component,  and  its  error  will  be  expressed  in  terms  of  an  expected  value; 
for  example,  the  RMS  error  in  latitude.  Note  that  if  there  are  no  parameter  errors,  the 
average  error  of  the  latitude  estimate  over  many  repetitions  of  the  experiment  is  zero7 
since  it  will  randomly  fluctuate  in  the  positive  and  negative  direction.  Of  course,  the 
mean  square  error,  which  is  the  expectation  of  a  positive  quantity,  is  still  nonzero.  Also, 
just  as  in  estimating  a  noisy  variable  by  averaging  over  a  set  of  samples,  for  a  single 
measurement  type  the  standard  deviation  of  the  state  estimate  will  vary  inversely  with 
the  square  root  of  the  number  of  measurements.  When  the  state  estimate  contains  more 
than  one  variable,  the  sensitivity  outputs  a  covariance  matrix,  i.e.,  E(AxsAxr),  where  ax8 
is  the  error  in  the  estimate  to  the  slh  variable.  For  example,  if  Xj  and  x2  represent  the 
latitude  and  longitude,  respectively,  then  we  obtain  E(Alat  Along)  as  well  as  E(Alat)2  and 
E(Along)2.  When  both  deterministic  and  stochastic  errors  are  present,  the  total  expected 
(RMS)  error  of  a  state  variable  is  computed  by  summing  the  deterministic  errors  and 
then  combining  that  result  with  the  other  errors  by  taking  the  square  root  of  the  surn  of 
their  squares. 

Figure  1.7  contains  an  example  of  the  output  of  the  software  package  SENS  for  a 
fixed  transmitter  and  a  moving  receiver.  For  completeness,  we  have  also  included  a  list¬ 
ing  of  the  SENS  input  requests  in  Figure  1.8.  In  the  case  illustrated,  the  unknown  state 
is  the  transmitter  latitude  and  longitude,  indicated  by  the  state  ty;,.;  STATEISLOC. 
The  orbit  model  is  the  NORAD  deep  space  model;  a  Kepler  model  has  also  been 

6  Stochastic  parameters  differ  from  measurements  only  in  their  functional  role  in  the  cost  function. 

They  enter  implicitly,  as  parameters,  and  thus,  conceptually,  are  not  as  central  to  the  estimation 
criterion. 

7  This  is  rigorously  true  only  to  the  first  order  since  our  nonlinear  estimates  are  •■nly  asymptotically 
unbiased. 
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doppler  type  is:  RELDOPPLER 
angle  type  is:  NONE 
element  estimates  type  is:  NONE 
State  type  is:  STATEISLOC 

Orbit  model  is:  NORAD 

Here  is  the  covariance  matrix: 

component  1  is  lat,  component2  is  long 
7. 420718742300576E+00  6 . 942098880222244E+01 

6. 942098880222244E+01  7 . 753117053419069E+02 


RMS  lit  error-  2.724099620480238E+00 
RMS  long  error-  2.784441964455188E+01 


Here  are  the  orbit  parameter  sensitivities: 
Units  are  degrees,  3cm,  and  seconds  of  time 


semi-maj  axis:  error  lat/long 
eccentricity:  error  lat/long 
epoch (2-body  rslt): 

error  lat/long 
mean  anomaly:  error  lat/long 
inclination :  error  lat/long 
ascending  node: error  lat/long 
arg  perigee:  error  lat/long 


LATITUDE  (deg) 

2. 653789595519299E-02 
-2. 061625284638502E+04 

5 . 423958191187525E-04 
-1. 298190677558637E-01 
5.973604795742479E-01 
-2. 559585668129556E-01 
6.899827141345592E-02 


LONGITUDE  (deg) 

2. 871645812066904E-01 
-1 . 092716754410192E+05 

4 . 273828602526120E-03 
-1 . 022914309755793E+00 
-8 . 025665754087354E+00 
-1. 563732042288754E+00 
1 . 692882756185627E+00 


Sensitivity  to  freq  params  per  Hz: 

LATITUDE  (deg) 

sat  offset:  error  lat/long  -  -5. 720171350036819E+00 

receiv  offset:  error  lat/long  -  5. 720171702619600E+00 

transmit  freq:  error  lat/long  -  -5 . 720171587202191E+00 


LONGITUDE  (deg) 

-6 . 170486123683860E+01 
6 . 170486225655878E+01 
-6 . 170486150408782E+01 


Sensitivity  to  receiver  location: 

gives  error  for  an  error  of  1  degree  in  indicated  parameter 
rev  latitude  (1  deg): 

tr  lat  (deg)-  6.819453815607811E-02 

tr  long  (deg)-  -3 . 260144317799733E-01 

rev  longitude  (1  deg): 

tr  lat  (deg)-  2.402443807119458E-01 

tr  long  (deg)-  2.420232834061336E+00 


Sensitivity  to  receiver  motion: 

gives  state  error  for  an  error  of  1  km/hr  in  spd  or  of  1  degree  in  course 
great  circle  speed  (1  km/hr): 


tr  lat  (deg) 
tr  long  (deg)- 
course  (1  deg)  : 
tr  lat  (deg)  - 
tr  long  (deg)« 


1. 199627 80727 3553 E+00 
1. 285460873073998E+01 

2 . 334775376994807E-01 
2 . 556731597113552E+00 


Figure  1.7a.  Example  of  the  output  of  the  sensitivity  software  (SENS)  where  the 
estimated  state  is  the  transmitter  position. 
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Sensitivity  to  satelite  position  per  RMS  meter 
(coords  are  in  geocentric  inertial  system) 

LATITUDE  (deg) 

xcoords  error  lat/long  -  1. 809036298711036E-04 

ycoords  error  lat/long  -  3 . 182725136278286E-04 

zcoord:  error  lat/long  -  1. 611327330902054E-05 


LONGITUDE  (deg) 

1. 344364853049826E-03 
3 . 493877932361680E-03 
1 . 124117784462831E-04 


Sensitivity  to  satelite  velocity  per  RMS  meter/sec 
(coords  are  in  geocentric  INERTIAL  system) * 

LATITUDE  (deg) 

xcoord:  error  lat/long  -  4.350484472181783E+00 

ycoord:  error  lat/long  -  2.490454480418472E+00 

zcoord:  error  lat/long  -  3.487122438460089E-01 


LONGITUDE  (deg) 

4 . 775662164106116E+01 
1. 850237658893264E+01 
3. 858319826910118E+00 


Figure  1.7  continued. 
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SENSITIVITY  INPUT  :  datafile 


give  element  epoch  (yr,day,  fractday) 

give  major  axis  {Jan}/  eccentricity,  meananom  (deg) 

give  inclination,  ascend  node,  argum  of  perigee  (deg) 

give  sn20  and  sn60  (revs/powerofday) 

give  bstar  (drag  term) 

give  iexp,  ibexp,  exponents  base  10  for  sn60  and  bstar 
give  offset  freqency  (MHz)  and  mixing  sign  (+-  1) 

give  station  epoch  (yr,day,  fractday) 
give  station  base  freq  (MHz) 

(this  is  transm  freq  or  receiver  offset) 
give  geodetic  latitude  and  longitude (deg) 
give  altitude  above  sea  (tan) 

give  motion  types  0  -  fixed/  1  -  greatcircle  motion 
give  course  (deg)  and  speed  (tan/hr) 
give  minimum  and  max  elevation  (deg) 

for  the  receiver  input: 
give  station  epoch  (yr,day,  fractday) 
give  station  base  freq  (MHz) 

(this  is  transm  freq  or  receiver  offset) 
give  geodetic  latitude  and  longitude [deg} 
give  altitude  above  sea  (tan) 

give  motion  type:  0  -  fixed;  1  -  greatcircle  motion 

give  course  (deg)  and  speed  (tan/hr) 

give  minimum  and  max  elevation  (deg) 

give  standard  dev  of  receiver  output  (Hz:  note  unit) 

for  doppler  measurements  give  (yr  day  fractday) : 
start  time 

upper  bound  (eg  stop)  time  for  meas 
give  interval  in  hrs  between  mesurements  and  max  no  of  pts 
for  azel  measurements  give  (yr  day  fractday): 
start  time 

upper  bound  (eg  stop)  time  for  meas 
give  interval  in  hrs  between  mesurements  and  max  no  of  pts 


SATELLITE 


TRANSMITTER 

if  motiontype  -  1 


RECEIVER 


if  motiontype  -  1 


MEAS  GENERATION 


Figure  1.8.  Input  requests  for  the  sensitivity  software  (SENS). 
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implemented,  and  other  models  may  be  interfaced  with  a  minimum  of  programming. 
The  type  RELDOPPLER  indicates  that  the  received  frequency  is  included  among  the 
input  measurements  in  the  form  of  equation  (1.16),  and  it  is  seen  that  its  standard  devi¬ 
ation  is  set  to  1.0  Hz.  Angle  type  NONE  implies  that  there  are  no  satellite  azimuth  or 
elevation  measurements.  The  satellite  orbit  information  may  be  in  the  form  of  parame¬ 
ters  (an  element  set)  or  ’’stochastic  parameters”  (position  and  velocity  at  each  time).  In 
other  words,  the  estimation  problem  corresponds  to  Figure  1.4. 

Q 

Each  of  the  entries  indicates  the  sensitivity  for  a  single  error.  To  evaluate  perfor¬ 
mance  when  there  is  more  than  one,  these  errors  must  be  combined  in  the  appropriate 
manner  as  described  at  the  end  of  the  previous  subsection.  Note  that  when  only  a  single 
type  of  input  error  is  present,  multiplying  that  error  by  a  factor  will  multiply  the  state 
estimate  error  by  the  same  factor;  however,  in  summing  different  types  of  errors,  the  sto¬ 
chastic  errors  must  be  combined  in  an  RMS  sense. 

The  first  group  of  results  shows  the  covariance  of  the  latitude,  longitude  for  a  1.0- 
Hz  standard  deviation  in  the  measured  receiver  frequencies.  We  see,  for  example,  that 
the  standard  deviation  of  the  estimate  of  the  latitude  is  2.7°.  Note  that  the  RMS  errors 
of  the  state  are  simply  the  square  roots  of  the  diagonal  elements  of  the  covariance 
matrix.  The  next  group  gives  the  parameter  sensitivities.  It  shows  the  errors  of  the  lati¬ 
tude  and  longitude  estimates  for  a  given  error  in  the  orbit  elements  or  in  the  frequency 
parameters.  The  input  errors  are  considered  to  be  a  single  unit  of  the  indicated  type. 
Thus,  we  see  displayed  respectively  the  lat/long  error  (not  RMS!)  for  a  1.0-km  error  in 
the  semimajor  axis;  for  a  1.0  error  in  eccentricity  (this  quantity  is  dimensionless);  for  a 
1.0-second  error  in  the  epoch  at  which  the  elements  are  supplied;  for  a  1.0°  error  in  the 
mean  anomaly,  inclination,  ascending  node,  or  argument  of  perigee;  and  for  a  1.0-Hz 
error  in  each  of  the  frequency  parameters.  The  next  two  groups  give  the  sensitivities  of 
the  state  to  errors  in  receiver  location  and  receiver  motion  respectively.  For  example,  an 
error  in  receiver  latitude  of  1.0°  results  in  an  error  in  the  estimated  transmitter  longitude 
of  -0.3°. 

The  last  group,  shown  in  Figure  1.7b,  gives  the  RMS  errors  in  latitude  and  longi¬ 
tude  as  a  function  of  errors  in  the  satellite  position  or  velocity.  In  other  words,  if 
instead  of  computing  the  satellite’s  orbit  via  an  orbit  model  driven  by  six  predetermined 
parameters  (elements),9  we  supply  a  set  of  satellite  location  and  velocity  ’’measurements” 
(assumed  to  be  corrupted  by  zero-mean  identically  distributed  noise),  then  the 
transmitter  location  will  have  the  indicated  RMS  errors.  This  is  particularly  useful  in 
determining  how  accurate  an  orbit  model  one  must  have  to  obtain  a  given  state  estimate 
accuracy.  (In  contrast,  the  element  sensitivities  cannot  take  into  account  any  errors  or 
approximations  which  may  reside  in  the  model  itself.)  The  output  is  arranged  analo¬ 
gously  to  that  above.  Thus,  for  example,  an  RMS  error  of  1.0  meter/sec  in  the  velocity 
of  the  satellite’s  y-coordinate  (i.e.,  the  values  supplied  for  vy  are  noisy  and  that  noise  is 
zero-mean  with  a  standard  deviation  of  1.0  meter/sec)  will  result  in  an  RMS  error  in 
latitude  of  2.5°  and  longitude  of  18.5°.  Actually,  the  entire  covariance  matrix  is 

8  An  exception  in  the  current  implementation  is  the  case  in  which  both  doppler  and 
azimuth/elevation  information  are  present.  In  that  situation,  the  covariance  matrix  combines  both 
(actually  three)  types  of  errors. 

*  In  the  NORAD  model  there  are  several  parameters  in  addition  to  the  six  elements,  but  we  do  not 
currently  compute  their  sensitivities. 
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computed  but  not  displayed. 

SUMMARY  OF  EQUATIONS 

We  make  the  following  definitions: 

x  =  N-dimensional  vector  whose  slh  component  is  the  sth  state  variable. 

x  =4  true  state  vector. 

xt  =  estimated  state  vector. 

ax  =  xt  —  x  the  error  in  the  state  estimate. 

q  =  an  arbitrary  parameter.  (That  is,  any  component  of  q.  The  remaining  components 

are  suppressed  for  notational  convenience.) 

A 

Aq  =  error  in  q. 

m(x,  q)  =  M-dimensional  vector  function  whose  ith  component  is  the  ith  measurement 
when  there  is  no  noise. 

m  =  vector  whose  components  are  the  actual  measured  values. 

We  note  that  the  cost  function  for  M  measurements  takes  the  form 


C(x)=E^4^, 

i=i 


and  the  estimate  xe  is  determined  by 


C(x,)  =  min  C(x). 


We  define  the  N  by  N  information  matrix  A  by 

M  3m:  i  3m: 

Ar,=  S' 


the  parameter  vector  b  by 


i=i  3xr  crj2  3xs  ’ 


m  3m:  1  3m: 

br  =  -  E' 


i=i  3xr  (Tj2  3q  ’ 


and  the  parameter  matrix  B  by 


„  ^1  ^mi  (  ^mi 

Brr' =  SiW 


(1.18) 


(1.19) 


(1.20) 


(1.21) 


(1.22) 
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A  numerical  procedure  for  the  solution  of  (1.18)  consists  of  the  iterations  (over  an 
index  k) 


xk+i  __  xk  _  ^  A-1VxC(xk) ,  (1.22b) 

where,  when  feasible,  X  =  1  (cf.,  Section  2). 

The  covariance  of  the  estimate  error  due  to  measurement  errors  (alone)  of  standard 
deviation  ax  is  computed  by 


X  «  A-1  . 


(1.23) 


Also,  for  a  single  parameter  error  of  Aq,  we  have  the  corresponding  state  estimate  error 

ax  *=»  A-1b  Aq,  (1.24) 

and  for  errors  in  satellite  position  or  velocity  (stochastic  parameters)  with  standard  devi¬ 
ation  crq,  we  have 


X  «  A^BA-Vq2  .  (1.25) 

Observe  that  the  matrix  A  is  a  sum  of  M  similar  terms  (equation  (1.20)).  Thus,  as 
the  number  of  measurements  increases,  holding  the  geometry  approximately  constant, 
the  entries  of  X  in  equation  (1.23)  will  vary  inversely  as  M.  We  may  interpret  this  as  a 
rule  of  thumb:  excluding  parameter  errors,  for  a  single  measurement  type,  the  standard 
deviations  of  the  errors  in  the  state  estimate  will  vary  inversely  as  the  square  root  of  the 
number  of  measurements  (l/\/M).  Note,  also,  that  the  matrix  in  equation  (1.25)  reduces 
to  products  of  the  components  of  (1.24)  if  we  substitute  Aq  for  crq.  The  approximation 
signs  are  used  rather  than  equalities  because  these  equations  are  only  correct  to  the 
first  order.  They  are  derived  in  Sect.on  2  by  linearization;  hence,  the  use  of  the  term 
sensitivity. 
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2.  SENSITIVITY  AND  ERROR  ANALYSIS 


DEFINITIONS 

In  the  previous  section,  we  discussed  a  nonlinear  least  squares  estimator  with  a  cost 
function  of  the  following  form: 


ow-  E 


i=l 


(2.1) 


where  x  is  the  state  vector,  rnj(x)  is  the  ith  measurement  function,  m;  is  the  ith  measure' 
ment,  and  <j\  is  its  variance.  The  state  estimate  xt  was  determined  by 


C(xe)  =  min  C(x)  .  (2.2) 

X 

Expression  (2.1)  is  appropriate  for  uncorrelated  measurements;  if  the  noise  also  happens 
to  be  Gaussian,  then  xe  is  the  maximum  likelihood  estimate  (MLE).  However,  since  gen¬ 
eralizing  to  correlated  measurements  does  not  complicate  our  derivations,  we  consider 
the  cost  function 


C(x)  =  (Am)*  R  -1(Am)  , 


(2.3) 


where 

f  =  matrix  transpose, 
x  =  N-dimensional  vector  of  state  variables. 

m  ==  (mj(x),  m^x),  •  •  •  ,  mM(x))*  column  vector  of  measurements 
m  =  vector  of  noisy  measurements 

=  m(x)  +  noise  where  x  is  the  true  state 
Am  =  m(x)  -  m 
R  =  correlation  matrix  of  Am. 

Note  that  in  the  uncorrelated  case,  the  components  of  m  are  the  m;  of  equation  (2.1), 
the  correlation  matrix  R  is  diagonal  with  entries  crj2,  and  (2.3)  reduces  to  (2.1). 

The  above  definitions  carry  the  implicit  restriction  on  the  cost  function  that  R  be 
the  correlation  matrix  of  Am.  This  insures  two  properties:  (i)  that  the  cost  function  be 
positive  semidefinite  and  (ii)  that  the  expected  value  E(im*  Am)  =  R.  The  second  con¬ 
dition  will  be  used  in  our  derivation  of  an  expression  for  the  sensitivity.  It  implies,  for 
example,  that  the  a  priori  variances  in  expression  (2.1)  be  chosen  with  reasonable 
accuracy.  The  resulting  simplification  is  well  worth  this  restriction.  For  most  applica¬ 
tions,  this  requirement  is  not  a  serious  problem.  The  fundamental  formula  (2.16),  which 
we  shall  derive,  may  be  modified  by  scaling  to  correct  approximately  for  a  least  squares 
weighting  which  is  not  exactly  the  inverse  of  the  covariance.  Note  that  if  <T\  is  a 
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constant  (i.e.,  if  the  measurement  noise  is  independent  and  stationary),  then  the  weight¬ 
ing  is  irrelevant  to  the  solution  of  (2.2). 

A  necessary  condition  for  (2.2)  to  be  satisfied  is  that  the  gradient  of  (2.3)  be  zero, 
i.e.,  that  xe  be  a  solution  of  the  set  of  N  nonlinear  equations  VC(xe)=0, 

Amt  R-i  4HL  =  0;  r  =  1, ...,  N.  (2.4) 


(We  have  used  the  symmetry  of  R  and  — —  =  — —  in  deriving  (2.4)  from  (2.3).) 

3xr  axr 

ERROR  SENSITIVITY 

We  wish  to  determine  the  accuracy  of  the  solution  xe  to  (2.4)  in  terms  of  the  accu¬ 
racy  of  the  measurements  m;  (the  ith  component  of  m)  and  the  parameters  which  enter 
into  the  function  m(-).  Let  us  first  consider  just  the  effect  of  measurement  errors.  To 
do  so  we  write 

m;  =  ihj  +  n; 

=  mi(x)  +  ni  ,  (2.5) 

where  n-,  is  the  the  noise  of  the  ith  measurement,  and  m;  is  the  value  of  the  measurement 
when  there  is  no  noise.  We  also  assume  that  the  measurements  are  unbiased;  i.e., 
E(n;)  =  0.  Define  a  vector  F  whose  components  Fr  are  the  left-hand  side  of  (2.4), 


Fr(xe,  m)  =  Am(xe,  m)f  R"1  -|^-(xe);  r  =  L  •••»  N. 


Then,  linearizing  about  (x,  m),  we  have 


N  dF 

Fr(xe,  in)  »  Fr(x,  m)  +  £  —  (x,  m)  (xe  -  x)4 

s=l  5xs 

M  <9Fr 

+  E  -to— (*» A)  ni  • 

i=l  ^mi 


Also,  since  Am(x,  m)  =  0, 


F r(x,  m)  =  0  , 

gFr  ,A  \ 9Amt  dm_ 

3xs  dx  8  oxT 


x=£ 


dFV 

dun 


ERii1 


3mj 

dxr 


(2.6) 


(2.7) 


(2.8a) 

(2.8b) 


(2.8c) 
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Thus,  to  the  first  order,  (2.4)  becomes 


E  R-1  f^x,  -  x).  =  n'  R-‘|!2-,  r  =  1 . N, 

s=i  oxs  3xr  3xr 


(2.9) 


where  the  partial  derivatives  are  understood  to  be  evaluated  at  x. 
For  convenience,  define  the  nxn  symmetric  matrix  A  by 


A  A  dm*  p_i  3m 
sr“  3xs  3xr 


(2.10) 


Letting  r  assume  two  values,  r  and  r' ,  and  multiplying  equation  (2.9)  by  itself,  we  obtain 


(A  (xe  -  x))r  (A  (xe  -  x))r.  =  Rij1  W  ■ 


(2.11) 


We  then  take  the  expectation  of  (2.11),  noting  that  E(njnji)  =  Rjj/  .  The  right-hand  side 
of  (2.11)  becomes 


e(rhs)  =  eAir,i,Al 


=  Arr»  ,  (2.12) 

where  E  denotes  expectation. 

The  expression  (xe  -  x)r  is  simply  the  rlh  component  of  the  error  in  the  estimate  xe. 
Thus,  X,  the  covariance  matrix  of  the  error  is 

X„  A  E  (x.  -  x)t(x,  -  i), ,  (2.13) 

and  substituting  (2.12)  into  the  expectation  of  (2.11),  we  have 

£Ars  ArV  Xss<  =  Arr!  ,  (2.14) 

ss' 

or 

AXA  =  A  (2.15) 

since  A  is  symmetric.  Solving  for  X,  we  find 
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X  =  A-!  . 


(2.16) 


Equation  (2.16)  along  with  definition  (2.10)  gives  a  simple  formula  for  the  covari¬ 
ance  of  the  estimate  xe.  (Actually,  since  it  was  derived  for  the  expected  squared  errors  it 
is  a  covariance  only  to  the  extent  that  the  estimate  is  unbiased.  It  is  easy  to  see,  by  tak¬ 
ing  the  expectation  of  (2.9),  that  xe  is  unbiased  up  to  the  first  order.)  If  the  noise  is 
Gaussian,  then  the  expression  (2.16)  represents  the  Cramer-Rao  bound  [3).  We  also  note 
that  the  diagonal  components  of  A-1  are  the  variances  of  the  components  of  the  state 
estimates.  In  practical  cases,  where  we  may  not  have  exact  knowledge  of  the  noise 
covariance,  it  is  often  a  reasonable  approximation  to  assume  that  we  know  it  up  to  some 
unknown  positive  factor  a-2,  i.e.,  E(n;nj)  =  R'jj/a2,  so  that  R  =  R'/a2.  Such  a  factor 
does  not  affect  the  solution  xe,  and  (2.16)  becomes 

X  =  a2A"1  .  (2.17) 

The  expected  cost  is  then  E(  C(x)  )  =  trace(R'  R-1)  =  Ma2;  thus,  we  may  approximate 
a 2  by  examining  the  average  of  the  residuals.  More  precisely,  the  sample  expectation  for 
C  is  given  by 


O  -  (m(xe)-m)t  R'1  (m(xe)-m)  , 


(2.18a) 


where  M  -  N  is  the  number  of  degrees  of  freedom  in  the  Gaussian  case.  (For  Gaussian 
noise,  C(x)  is  a  x2  distribution.)  We  then  estimate  a  by 

a2«C/M.  (2.18b) 


Finally,  we  note  that  if  the  measurements  are  independent  (i.e.,  Rjj  =  /  o\  ), 

then  (2.10)  becomes 


or 


M  dm;  am;  1 
saxr  sx3  ^i2 

(2.19) 

^  Vm-,  ^  Vmj 
=  h  ®  > 

i=i  <*i  <7i 

(2.20) 

where  @  stands  for  the  matrix  tensor  product. 

Next  we  compute  the  sensitivity  to  deterministic  parameter  errors.  The  procedure 
is  almost  identical  to  that  which  led  to  equation  (2.16).  Let  q  stand  for  some  parameter, 
so  that  the  measurement  functions  mj(x,  q)  are  dependent  on  the  state  x  and  a  parame¬ 
ter  q.  Previously  the  dependence  on  q  was  left  implicit.  Similarly,  we  now  consider  only 
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one  parameter  q  at  a  time  and  set  all  other  errors  equal  to  zero.  (Since  the  sensitivity  is 
a  linearization,  the  first-order  error  in  x  due  to  all  the  errors  will  be  the  sum  of  the  indi¬ 
vidual  contributions.)  In  the  absence  of  measurement  noise  we  write  (cf.,  equation  (2.5)) 

mj  =  mi(x,  q)  ,  (2.21) 

where  q  is  the  true  value  of  q.  As  in  (2.6)  we  define 

Fr(x*>  q)  =  Am(xe,  q)t  R"1  -|^-(xe,  q).  (2.22) 

Then,  linearizing  about  (x,  q),  we  have 

N  3F, 

Fr(*e,  q)  »  Fr(x,  q)  +  E  q)  (x*  -  *)* 

5=1  °XS 

+  q)  (q  -  q)  •  (2-23) 

dq 


Also,  since  Am(x,  q)  =  0, 


aFr 

dx. 


Fr(*,  q)  =  0 
..  dAm1  D_,  9m 

(x' q)  “  UT  R  &T" 

am 

dxr 


5Fr  A  aAm^  -p-i 

-ar(x' q) = R 


q=q 

x=x 


q=q 

x=x 


Thus,  to  the  first  order,  (2.4)  becomes 


or 

A  (xe  -  x)  =  b  Aq  , 

where 

■  A  am*  tj-1  dm 

h' - 9TR  "9x7 


(2.24aj 

(2.24b) 

(2.24c) 


r  =  1,  •••,  N 


(2.25) 


(2.26a) 

(2.26b) 


-  q  • 


(2.26c) 


25 


Note  that  if  R  is  diagonal,  then  (2.26b)  reduces  to  (cl\,  equation  (2.19)) 


m  am.  am.  j 
kf  Zj' 


i»i  dq  dxT  o? 


Solving  (2.26a)  for  ax  =  xe  -  x,  the  error  in  x,  we  get 

ax  =  A_1b  Aq  . 


(2.26d) 


(2.27) 


Equation  (2.27)  is  an  expression  for  the  error  (to  the  first  order)  in  the  estimate  xe 
caused  by  a  single  parameter  error  Aq,  whereas  (2.16)  gives  the  mean  squared  error 
E(Axr)2  =  Aji1  due  to  an  entire  set  of  noisy  measurements.  The  derivation  of  (2.27) 
does  not  require  that  the  least  square  weights  R"1  be  the  actual  covariance  of  the  meas¬ 
urements  rh;.  In  contrast,  expression  (2.16)  was  derived  under  the  assumption  that 
E(n;nj)  =  Rjj. 

Finally,  let  us  consider  the  case  in  which  some  of  the  parameters  are  themselves 
random  variables,  i.e.,  where  the  measurement  function  is  of  the  form  m;(x,  qa),  where 

qa  =  qa+na  ;  E(na)  =  0  ;  a  =  1,  2,  •  •  •  (2.28) 

are  a  set  of  noisy  parameter  estimates.  The  prototype  for  this  in  the  current  study  is  a 
set  of  satellite  positions  and  velocities.  Note  that  ir  these  are  supplied  for  each  measure¬ 
ment  time  m„  then  the  index  a  will  run  from  1  to  6M,  and  m;(x,  qa)  will  only  be  a  func¬ 
tion  of  qa  for  the  six  parameters  relevant  to  that  time,  e.g.,  a  =  6(i  -  1)  +  1  to  a  =  6i. 
This  property  implies 


3m;  dmj» 

3qa  dqa 


if  i  7^  i' . 


(2.29) 


Observe  that  the  only  difference  between  the  ’’stochastic  parameters”  of  equation  (2.28) 
and  the  noisy  measurements  of  equation  (2.5)  is  the  functional  manner  in  which  they 
enter  the  cost  function. 

Continuing  with  our  development,  we  note  that  equation  (2.26)  remains  valid 
except  that  q  is  replaced  with  a  sum  over  qa  . 


(AK-x))r=-S^R^ 


Aqa  ;  r  =  1,  ...,  N 


aij 


(2.30) 


As  in  the  derivation  of  (2.14),  we  take  two  values  of  r,  multiply  (2.30)  by  itself,  and  take 
expectations.  If  we  assume  that  the  noisy  parameters  are  independent  and  identically 
distributed, 
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we  get 


(2.31) 


E(AqaAqa, )  =  <rq2  Sia!  , 


S^rs^-rV  Xss' 
ss' 


_  2v-i  ^ dmi 
Zj  L  on 

a  iji'j' 


3nij 

3xr 


3m ;» 
3q& 


3my 

3xr/ 


This  may  be  written  as 


AXAf  =  B  <Tq2  , 


(2.32) 


(2.33) 


with  B  defined  by  (2.32).  If  R  is  diagonal,  B  simplifies  to 


Brr'  =SS 


3nij  i  3m;  3m;i  i  dm;/ 


ii'  5xr  o\  d<h  5cia  o?  3xri 


Furthermore,  if  equation  (2.29)  holds,  then 


„  3m: 

Brr'  =  S  S— T  "fa- 

a  i  ff\  oXf 


dm\  \  dm\ 


3qa  dxr» 


In  any  case,  the  solution  to  (2.33)  is  given  by 


X  =  A-1BA-W 


(2.34) 


(2.35) 


(2.36) 


since  all  the  matrices  involved  are  symmetric. 

Inasmuch  as  the  above  computations  are  first-order  approximations 
(and  E(n;)  =  0),  it  is  easy  to  see  that  if  we  assume  the  various  stochastic  types  of  errors 
are  independent,  then 

(1)  Parameter  errors  add  when  more  than  one  is  present. 

(2)  The  total  stochastic  error  is  the  RMS  error  of  the  individual  errors. 

(3)  The  total  error  is  the  RMS  sum  of  the  errors  in  (1)  and  (2). 

SOLVING  FOR  THE  STATE:  COMPUTATIONAL  CONSIDERATIONS 

The  major  computational  steps  in  the  above  sensitivity  analysis  involve  finding  (a) 
the  gradients  of  the  measurement  function,  3mj/3xr  and  3mj/3q,  (b)  the  matrix  A  which 
equals  their  outer  product,  and  (c)  the  inverse  of  A.  Typically,  (b)  and  (c)  are  inexpen¬ 
sive  in  comparison  with  (a).  In  light  of  this,  and  noting  that  the  most  common  methods 
for  solving  (2.2)  are  steepest  descent  methods  utilizing  the  gradient  of  the  cost  function 
(VC  =  Am’  R"1  Vm),  we  see  that  the  same  computations  are  involved  in  determining 
the  estimate  xe  and  the  sensitivity.  In  other  words,  the  sensitivity  can  be  computed  at 
essentially  no  additional  cost. 
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Finally,  we  note  that  a  particularly  effective  steepest  descent  algorithm  is  the 
modified  Newton-Raphson  algorithm 


x(i+l)  =  x(i)  -  A-1(i)  VC(i) 


or  the  more  general 


x(i+l)  =  x(i)  -  X(i)  A-1(i)  VC(i)  . 


(2.37) 


(2.38) 


Since  A  is  positive  definite  unless  the  state  x  is  unobservable  (i.e.,  unless  the  system  of 
equations  (2.4)  is  singular),  the  direction  A-1  VC  is  one  of  decreasing  cost,  and  (2.38) 
will  converge  for  suitable  choices  of  the  step  size  X(i).  Recursion  (2.37)  is  simply  an  N- 
dimensional  Newton-Raphson  iteration  for  solving  the  nonlinear  system  of  equations 
(2.4)  but  with  the  Jacobian  of  that  system, 


J 


rs  — 


dAm*  p_i  dAm 
dxr  dxs 


+  Am*  R"1 


d2m 

dxrdxs  ’ 


(2.39) 


replaced  by 


dAm*  dAm 

dxt  dxs 

That  is,  we  use  the  matrix  A  in  place  of  the  Jacobian.  Such  an  approximation  is  intui¬ 
tively  reasonable  since  for  x  =  x  with  zero  noise,  Am  =  0  and  the  Jacobian  actually 
equals  A.  Furthermore,  even  if  it  happens  that  the  cost  does  not  decrease  during  an 
iteration  of  the  form  (2.37),  it  is  still  guaranteed  to  do  so  (unless  one  has  arrived  at  the 
minimum  xe)  for  some  X  <  1  in  (2.38). 

COMPOSITION  OF  EXPERIMENTS 

Suppose  that  two  experiments  are  performed  in  succession.  The  first  is  used  to 
determine  a  set  of  satellite  orbit  elements  ys,  where  s  =  1,  ...,  Nj  (e.g.,  Nj  =  6  ).  Let  its 
covariance  matrix  be  represented  by  Y.  The  second  experiment  employs  these  elements 
as  parameters  in  performing  a  localization,  that  is,  to  determine  a  position  xr  where  r  = 
1,  ...,  N2  (e.g.,  N2  =  2).  We  wish  to  determine  the  covariance  of  x  due  to  the  orbit 
errors  of  the  first  experiment.  From  (2.27),  the  errors  in  x  are  (deterministically)  related 
to  each  of  those  in  y  by 


Ax(due  to  ys)  =  A-1bs  Ays  , 


(2.41) 


where  bs  is  defined  to  be  b  of  equation  (2.26b)  with  q  replaced  by  ys.  If  we  let 
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D„  4  (A-ib*)t 


(2.42) 


be  the  matrix  whose  columns  are  the  sensitivity  of  x  to  yS)  then  the  sensitivity  to  all  the 
errors  is  given  by  the  their  sum,  i.e.,  by  the  sum  over  s  of  the  right-hand  side  of  (2.41). 
This  yields 


Ax  =  DAy  .  (2-43) 

Finally,  we  have  for  the  covariance  of  x  due  only  to  the  errors  in  y 

E(AXAxt)  =  E(DAy(Ay)tD‘) 


=  DYDf  .  (2.44) 

It  should  be  emphasized  that  if  the  two  experiments  are  not  exactly  the  same,  the  meas¬ 
urement  functions  used  in  determining  Y  and  those  in  D  may  not  be  the  same. 

Let  us  examine  (2.44)  a  little  more  closely.  Suppose  that  the  measurement  spaces  of 
the  two  experiments  coincide.  Define  a  set  of  Nj  vectors  ks  and  a  set  of  N2  vectors  xr  in 
Euclidian  M-dimensional  space  by 


a  9m  ,  NI  r  3m  ,  KT 

* “  If,  s  =  1 . N,i  x  a?  r  =  1 . Ns- 


(2.45) 


Thus,  for  example,  the  ith  component  of  #c*  is  equal  to  dm\/dy&.  We  define  an  inner  pro¬ 
duct  on  that  space  by 


Then  we  have 


m  m'  =  Yj  niiRi^m'j  . 
U=i 


M 


(2.46) 


Arr'  =  Xr  ■  XT> 

(2.47a) 

(bS)r  =  -Xf  • 

(2.47b) 

^  -  Xs  •  k*'  . 

(2.47c) 

The  matrix  D  becomes 
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(2.48) 


Drs  =  EAfk1  (-Xk-«s), 
k 

and 

(DYD'y  =  EA,i;1(xkOYm„(*"-Xi)AiTf.  (2.49) 

kmni 

Let  P  be  the  operator  defined  on  an  arbitrary  vector  a  by 

P«=E«mY>‘4  (2-5°) 

mn 

Then  P  is  linear  and  k  projects  vectors  onto  the  space  spanned  by  the  k3,  s  ==  1,  Nj. 
To  see  this  we  note  that  P/c^  =  5jKmYmn(«n-K^)  =  ^]/cmYmnYn'j1  =  «:■*,  and  if  a-/cm  =  0 

mn  mn 

for  all  m,  then  Pa;  =  0.  To  take  advantage  of  the  operator  P  in  equation  (2.49),  we 
decompose  the  vectors  xr  into 


Xr  =  X^  +  X- 


(2.51) 


where  x+  is  contained  in  the  space  spanned  by  the  n'  s,  and  X-  is  orthogonal  to  it. 
Note  that  the  orthogonality  produces  a  corresponding  decomposition  in  A: 

V  =  xl-xt  +  Xl-xl' 

=  (A+V  +(A.)r,  .  (2.52) 

Equation  (2.49)  becomes 

(DYD')„.  =  EAS‘(xk-PxW 

ki 

=  EArkHx+'xijAiT'1  (2.53) 

ki 


or 

DYD*  =  A"1  (A  -  A.)  A"1 


=  A-1  -  A-1  A. A-1. 


(2.54) 
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Since  A  and  A_  are  positive  semidefinite10  and  since  A"1  =  (A'1)',  we  conclude  that 


DYD*  <  A"1 


(2.55) 


as  quadratic  forms.  In  other  words,  if  the  measurement  functions  are  the  same  in  the 
two  experiments,  then  the  influence  of  the  errors  in  the  first  experiment  on  the  covari¬ 
ance  of  the  output  of  the  second  is  at  most  as  great  as  the  influence  of  the  measurement 
errors  in  the  second  experiment  on  its  output: 

E(axAX*)  <  A"1 
The  total  error,  if  both  experiments  are  noisy,  satisfies 

Etot^AxAx1)  <  2A'1 .  (2.57) 


(2.56) 


10  To  see  this,  note  that  they  are  sums  of  outer  products  of  vectors.  For  example, 
£)z rAn/  —  S^X^X1'  Zf  =  |  w  |  2  where  w  =  2zrXr  ^or  an  arbitrary  Nydimensional  vev‘or  z. 
rf  r 

Hence,  A  is  a  nonnegative  quadratic  form. 
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3.  COMPUTATIONS  IN  GEOCENTRIC  COORDINATES 


In  order  to  implement  the  sensitivity  equations  of  Section  2,  as  well  as  to  solve  the 
estimation  problem  via  (2.38),  we  require  the  gradients  of  the  measurement  functions 
fft(x,  q,  t)  or  d(x,  q,  t)  (cf.  equations  (1.13)  and  (1.14))  with  respect  to  the  variables  x 
and  q.  Since  distinguishing  a  variable  as  being  part  of  the  state  variable  x  or  one  of  the 
a  priori  parameters  q  is  irrelevant  to  such  a  computation,  we  simply  employ  the  notation 
q  for  the  independent  variable  in  the  derivative.  On  the  other  hand,  the  difference 
between  a  ground  station  parameter  and  a  satellite  parameter  will  prove  useful,  and  for 
these  we  shall  write  qc  and  qs,  respectively. 

We  begin  Section  3  with  a  computation  of  the  derivatives  of  the  relative  doppler 
function  b  ( cf.,  equation  (1.8))  with  respect  to  qc  and  qs.  This  brief  derivation  is  at  an 
intermediate  level  of  detail,  almost  independent  of  the  coordinate  system.  In  the  next 
subsection  we  introduce  standard  geophysical  coordinates  and  describe  the  computations 
leading  from  latitude,  longitude,  and  velocity  on  an  oblate  earth  to  the  vectors  R,  r,  R, 
and  r  of  Figure  1.3.  We  also  provide  equations  for  satellite  azimuth  and  elevation  meas¬ 
urements.  The  extension  to  ground  station  motion  (i.e.,  great  circle  motion  parameter¬ 
ized  by  initial  course  and  speed)  is  given  in  Section  4. 

DOPPLER  DERIVATIVES 

We  proceed  to  compute  the  derivatives  of  the  relative  doppler  function  b  of  equa¬ 
tion  (1.8).  Noting  that 


dq  ||  p  ||  dq 


-  -sr(P  ’  p) 


,-1/2 


we  have 


d_ 

dq 


IIP  II 


dp 

dq 


P  ,  dp 
I!  P  II  dq 


d  •  dP- 

P  dq 
II  P  II  3 


P  (P-P) 

Up  II  II p  II ; 


(3.1) 


(3.2) 


Thus,  if  we  define  the  normal  vector  n  by 


n  — 


Up  II  ’ 


then 


db_  _  1  dp  _  1  dp 
dq  c  dq  n  ||  p  ||  dq 


A_  (P'PLn 

c  c  II  p  II 


(3.3) 


1  dp  _  1  dp 

c  dq  ”  ||  p  ||  dq 


A  +  b  n 

c 


(3.4) 
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As  a  consequence  of  (1.6)  and  (1.9),  equation  (3.4)  is  a  function  of  r,  r,  R,  and  R. 
We  can  remove  the  dependency  on  R  by  taking  into  account  the  earth’s  motion.  Let  us 
denote  by  Cl  the  vector 


Cl  =  Wq  K 


(3.5) 


where  K  is  a  unit  vector  pointing  north  along  the  earth’s  axis,  and  w0  is  the  rate  of  the 
earth’s  rotation  (e.g.,  radians/sec).  Then,  if  a  ground  station  is  motionless  relative  to 
the  earth  with  position  vector  R  (cf.  Figure  1.3),  its  velocity  with  respect  to  any  geocen¬ 
tric  inertial  coordinate  syste7n  (cf.  next  subsection)  is  R  =  Cl  x  R,  where  x  indicates  the 
vector  product.  Finally,  if  the  station  is  moving  relative  to  the  earth  with  a  ’’ground 
velocity  of  v,  we  have 


R  =  Cl  x  R  4-  v  . 


(3.6) 


That  is,  v  denotes  the  velocity  the  station  would  have  if  we  instantaneously  stopped  the 
earth’s  rotation.  Equation  (1.9)  becomes 


p  =  r  -  O  x  R  -  v  . 


(3.7) 


The  derivatives  of  b  simplify  if  we  consider  satellite  parameters  and  ground  station 

0R  dR 

parameters  separately.  For  the  satellite  ==  —  0  so  that 


•  n _ _ _ —  •  U  +  bn 

3qs  c  dqs  [|  P  ||  <9qs  c 


(3.8) 


For  the  ground  station 


-21—  =  -£L-  =  o  and  using  (3.6)  with  (3.7), 
3qG  dqG 


db 

5qc 


i(nx®  + 

c  c>qG  dqG  J 


1  8R 

IIP  II  3qG 


dR 

dqo 


n  x  Cl  ,  p 
c  c  ||  p  1] 


+ 


b  n 

IIP  II  . 


1  dv 
c  dqG 


n 


_1 _ dR 

Up II  5cig 


\pxJL  +  ±  +  bn 

{  c  c 


1  dv 
c  dqG 


n  . 


(3.9) 
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(We  have  used  the  vector  indentity  (A  x  B)  ■  C  =  B  •  (C  x  A). 

For  completeness,  we  compute  the  derivatives  of  the  measurement  function 
d(x,  q,  tj)  with  respect  to  the  various  parameters.  If  q  is  not  a  frequency  parameter  then 
we  have 


Also,  . 


ad 

aq 


abu  bD  +  ^(bu  +  1  +  ip.) 

In 


aq 


aq 


for  q  7^  freq  param  . 


ad 

dgB 


=  1 


y 


(3.10) 


(3.11) 


and 


where 


ad 


=  j-(bD  +  f) , 
!0 


(3.12) 


ad  £_ 

dfRO  fo 


(3.13) 


{0  for  gg  a  state  variable 
1  otherwise 

The  partial  derivative  with  respect  to  f0  which  is  suitable  for  use  in  the  sensitivity  is 
somewhat  more  complex,  since  the  modified  measurement  d;  of  (1.16)  is  a  function  of  f0. 
One  can  either  work  entirely  in  terms  of  the  measurement  functions  fR(x,  q,  t;)  of  (1.13) 
or  utilize  the  derivative  of  the  entire  residual  (which  appears  in  (1.17)), 


d(d  -  dj) 

af0 


==  -72-(bDfsO  +  e(fSO  -  fRo)  “  ?i) 

‘0 


(3.14) 


GEOCENTRIC  COORDINATES 

Equations  (1.6),  (1.8),  and  (3.7)-(3.14)  determine  the  measurements  at  an  arbitrary 
time  t  as  a  function  of  the  vectors  R(t),  r(t),  r(t),  and  v(t)  at  that  time.  However,  the 
parameters  q  are  fixed  numbers  associated,  at  least  conceptually,  with  some  fixed  time  t0. 
Even  for  a  stationary  ground  station,  it  is  natural  to  equate  the  latitude  and  longitude 
at  some  time  t0  with  a  fixed  vector  R0  and  then  to  compute  R(t)  as  a  function  of  Rq  and 
t  —  t0  by  taking  into  account  the  earth’s  rotation. 

In  order  to  be  more  precise,  let  us  introduce  some  specific  coordinate  systems.  Fig¬ 
ure  3.1  illustrates  the  geocentric,  or  what  we  shall  call  the  IJK  coordinate  system.  The 
K-axis  is  the  axis  of  the  earth’s  rotation;  the  other  two  axes  are  chosen  so  that  IJK 
forms  a  right-handed  system  of  mutually  perpendicular  axes  with  I  pointing  in  the 
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direction  of  the  vernal  equinox11  [l,  2).  The  vernal  equinox  and  the  earth’s  axis  both 
move  with  respect  to  the  celestial  sphere.  Thus,  to  specify  an  UK  coordinated  system 
we  must  also  supply  an  epoch  (point  in  time)  r.  (Note  that  t  is  yet  a  third  time,  that  of 
the  measurement.  It  usually  is,  but  need  not  be,  the  same  as  t.)  A  typical  computation 
involves  the  determination  of  the  satellite  position  r(t)  in  IJK  coordinates  relative  to  one 
epoch,  rs  and  the  location  of  a  ground  station  R(t)  in  a  coordinate  system  of  another 
epoch  tq.  Since  vector  computations  such  as  p  =  r(t)  -  R(t)  must  be  done  in  the  same 
coordinate  system,  we  must  either  transform  R  from  IJKG  to  IJKg,  or  r  from  IJKS  to 
IJKq,  or  both  vectors  to  some  other,  common  coordinate  system  before  subtracting. 
Such  a  transformation,  between  orthogonal  coordinate  systems,  will  be  a  rotation,  which 
may  be  identified  with  a  matrix  T.  We  use  the  notation 

a(s)  =  T(rs,  rG)a(G) 


=  TSGa(G>  (3.15) 


for  the  mapping  of  the  representation  of  an  arbitrary  vector  a  in  IJKG  coordinates  to 
IJKs  coordinates.  One  should  be  careful  not  to  confuse  vectors  as  abstract  entities  with 
their  representations;  thus,  a  is  a  single  object,  a(s)  is  its  representation  as  a  column  of 
three  numbers  in  the  IJKS  coordinate  system,  while  a^G^  is  a  representation  of  the  same 
object  in  the  IJKG  coordinate  system.  The  matrix  T  is  a  function  only  of  the  two  epochs 
rs  and  rG;  the  details  may  be  found  in  (2).  Finally,  we  have 

p(S)  =  r(S)_R(S) 


=  r<s>  -  TSGR(G)  (3.16) 

p(s)  =  f.(S)  _  Tsg  (fi(G)  x  R(G)  _  v(G)j  =  -  TsgR(g)  .  (3.17) 

The  computation  of  R(g)  is  detailed  below;  that  of  and  are  described  in  Appen¬ 
dix  B. 

From  Figure  3.2,  we  see  that  the  IJK  coordinates  of  a  station  are  determined  by  0, 
the  angle  of  its  meridian  with  respect  to  the  I-axis,  and  by  its  rectangular  coordinates 
(x,  z)  in  the  plane  of  that  meridian.  That  is, 


R(g)  =  (x  cos0,  x  sin0,  z)  (3.18) 


11  The  vernal  equinox  is  the  intersection  of  the  equatorial  plane  (the  IJ  plane  through  the  geo¬ 
center)  with  the  ecliptic. 
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K 


Figure  3.2.  A  meridial  slice  through  a  ground  station  at  R.  For  an 
oblate  earth  this  slice  is  an  ellipse.  See  Figures  3.3  and  3.4. 


According  to  Figure  3.1,  6  is  simply 


*00  =  0,(t)  +  X  , 


(3.10) 


where  X  is  the  longitude  of  the  ground  station,  0g(t)  is  the  angle  of  the  Greenwich  meri¬ 
dian  with  respect  to  I  at  time  t  (i.e,,  the  greenwich  sidereal  time),  and  coordinate  system 
IJKq  is  that  of  the  epoch  rG  =  t. 

The  earth  is  more  accurately  represented  by  an  ellipsoid  than  the  sphere  drawn  in 
Figure  3.1.  In  that  case  the  meridial  slice  is  an  ellipse,  as  pictured  in  Figure  3.3.  The 
geodetic  latitude  shown  there  is  determined  by  the  angle  subtended  by  a  line  perpen¬ 
dicular  to  the  earth’s  surface.  Some  straightforward  algebra  yields  the  rectangular  coor¬ 
dinates  of  the  ground  station  [l], 


x  = 


t  Vl  -  e2sin 


+  H  cos  <f> 


-(  ^-e)  ■« 

\  \/l  -  e2sin2<£ 


smfi 


(3.20) 


where  the  ellipsoid  is  taken  to  lie  at  mean  sea  level,  H  is  the  altitude  (i.e.,  height  of  the 
ground  station  above  sea  level),  ae  is  the  equatorial  radius  of  the  earth,  and  e  is  the 
earth’s  eccentricity.  Equations  (3.18)  to  (3.20)  determine  R^(t),  given  the  latitude  and 
longitude  <j>  and  X  for  a  motionless  ground  station  (i.e.,  v  ==  0). 

The  computation  of  v(G)(t),  the  instantaneous  ground  station  velocity  relative  to 
the  earth  in  IJK  coordinates,  is  not  quite  so  straightforward.  If  we  are  given  a  course 
iv(t)  and  speed  s(t)  for  each  time  t,  then  the  procedure  is  a  simple  one,  which  we  describe 
below.  However,  if  we  wish  to  model  the  motion  parametrically,  for  example,  in  terms  of 
a  constant  speed  and  an  initial  course  at  some  time  t0,  then  we  must  specify  a  curve. 
We  cannot  simply  use  a  rectilinear  motion  of  the  form  a  =  a0  +  v0  (t  -  t0)  if  we  wish 
the  station  to  remain  on  the  surface  of  the  earth.  A  natural  choice  for  a  spherical  earth 
is  great  circle  motion,  i.e.,  motion  along  a  great  circle  with  constant  speed  and  an  initial 
course  vQ.  Note,  that  in  that  case,  J^t)  ^  v0.  Motion  on  an  oblate  earth  is  even  less 
obvious.  These  are  treated  in  Section  4. 


AZIMUTH  AND  ELEVATION 

In  addition  to  the  IJK  coordinates,  Figure  3.1  illustrates  the  so-called  SEZ  coordi¬ 
nate  system.  The  letters  S  and  E  stand  for  south  and  east.  The  Z-axis  is  a  perpendicu¬ 
lar  to  the  earth’s  surface  at  the  ground  station;  S  and  E  are  placed  in  the  plare  tangent 
to  the  earth.  (Note  that  these  definitions  are  valid  for  an  oblate  earth  even  though  Fig¬ 
ure  3.1  is  misleading.)  SEZ  coordinates  are  fairly  natural  for  describing  local  motion. 
For  example,  an  object  moving  at  a  speed  s  and  course  v  (clockwise  from  north)  has  a 
velocity  with  SEZ  coordinates 
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normal 


Figure  3.3.  Illustration  of  the  definition  of  geodetic  latitude. 
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-s  cosy 


V(SEZ)  = 


s  siny 
0 


(3.21) 


This  SEZ  coordinate  system,  of  course,  moves  with  time  as  the  earth  rotates.  It  also 
depends  on  the  ground  station  location,  i.e.,  on  <f>  and  X.  The  transformation  from  SEZ 
coordinates  at  time  t  to  IJK  coordinates  at  the  same  epoch  (i.e.,  r  =  t)  is  given  by  [l] 


where 


a(IJK)  = 

D-1  a(SEZ> 

> 

(3.22a) 

sin<£  cos# 

sin  if)  sin# 

-cos^ 

-sin# 

cos# 

0 

(3.22b) 

COS0  cos# 

cos^  sin# 

sinii 

(The  fact  that  this  transformation  holds  even  though  it  does  not  reflect  the  oblateness  of 
the  earth  is  clarified  in  Figure  3.4.  If  the  origin  of  the  IJK  system  is  translated  from  A 
to  B,  then  the  transformation  is  the  same  as  would  be  obtained  in  Figure  3.1.  Since  we 
are  describing  only  detached  vectors,  i.e.,  a  length  and  a  direction,  in  this  transformation 
the  translation  from  A  to  B  has  no  effect.) 

The  relationship  of  the  so-called  topocentric  coordinate  system  (i.e.,  azimuth  and 
elevation)  to  SEZ  coordinates  is  described  by  Figure  3.5.  From  that  figure, 


AZ  =  -  tan-1— 

Ps 

EL=si”"‘w  (3'23) 

with  the  sign  of  AZ  determined  by 

sin  AZ  =  ^=p=¥  ■  (3.21) 

V  Pe  +  Ps 

Next  we  compute  the  derivatives.  Using  (tan-1x)'  -'<=  1/(1  4-  x2),  we  have 

5AZ  =  Ps2  d  Pe 
5cl  Pe  +  Ps2  5cl  Ps 

Ps2  (  1  1  f  dpE  <9ps  ) 

~  Pe  +  Ps2  [  Ps2  )  (  PS  30  "  PE  30  J 
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1 

Pe  +  Ps2 


(3.25) 


An  alternative  form  obtained  by  using  n  =  p/  ||  p  ||  (definition  (3.3))  and  +  Pe  = 
||p||2(l-nz2)  is 


dAZ 

dq 


II  P  II  (l~nz2)  ‘ 


Also,  since 


sin  EL  =  nz  , 


dEL 

3q 


cos  EL  = 


dnz 

<9q 


Then,  cos  EL  =  y/l  -  nz  implies 


dEL 

dq 


dnz 

.  ga  . 
V1  ~  nZ 


Then,  since 


dnz 

dq 


3pz 


dq 

II  pH 


we  have 


1 

II  pH 


dpz  dp 
dq  nz  dq 


dEL 

dq 


dpz 

"aq”"z 


II  P  II  V1  -  nZ 


(3.26) 


(3.27) 

(3.28) 


(3.29) 


(3.30) 


(3.31) 
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TIME  AS  A  PARAMETER 

It  is  readily  apparent  from  the  current  discussion  that  timing  represents  a  potential 
physical  source  of  error  in  the  geolocation  computations.  Table  3.1  below  summarizes 
the  relevant  variables. 


Table  3.1.  Summary  of  time  variables. 


initial  state 

signal  time 

pos/vel 

coordinate  sys 

receiver 

tR0 

t 

RR(fc)>  RR(t) 

% 

transm 

tT0 

t 

R-r(t),  R-r(t) 

rT 

satellite 

tso 

t 

r(t),  r(t) 

rs 

The  parameters  tR0  and  tT0  represent  the  times  at  which  the  ground  station  states  are 
specified.  Similarly,  tso  stands  for  the  epoch  for  the  satellite  elements,  t  is  time  of  the 
receiver  measurement;  hence,  effectively,  it  also  equals  the  time  of  the  signal  transmission 
and  of  its  relay  at  the  satellite,  tq  and  r$  are  intermediate  variables,  not  parameters; 
that  is,  they  are  either  functions  of  parameters  or  errorless  constants.  In  the  present 
case,  rR  =  rT  =  t,  while  rs  is  the  epoch  which  the  satellite  model  uses  for  its  coordinate 
system.  Usually  rs  =  tso,  the  epoch  of  the  element  set  (Kepler)  or  some  fixed  astronomi¬ 
cal  date  as  in  ts  =  Besselian  year  1950  (NORAD). 

The  parameters  tR0  and  tTo  only  have  significance  when  a  ground  station  is  moving. 
Even  in  that  case,  an  error  in  t0  may  be  equated  to  an  appropriate  error  in  the  initial 
latitude,  longitude,  course,  and  speed  parameters.  Similarly,  tso  is  equivalent  to  a  set  of 
errors  in  the  satellite  parameters.  For  a  Kepler  model,  it  is  equivalent  to  an  error  in  the 
mean  anomaly.  Of  greater  interest,  perhaps,  are  errors  in  t,  the  measurement  times. 
These  may  usually  be  modeled  as  a  bias  in  the  setting  of  the  receiver’s  clock.  More 
specifically,  we  write 


t  —  t*  +  CR  , 


(3.32) 


where  cR  is  the  clock  offset.  Then  an  error  in  cR  corresponds  to  a  bias  in  the  clock.  The 
sensitivity  to  this  bias  is  a  function  of  the  partial  derivatives  with  respect  to  cR.  From 
(3.16),  (3.17),  and  (3.22),  we  have 


dp(s)  _  dr^  _  aTscR(G) _  ^  dR^ 

<9cr  dt  dt  SG  dt 


(3.33a) 


=  _  dTso  x  r(g)  _  V(G)| 

<9cr  dt  dt 


44 


(3.33b) 


-Tsg 


n<°>  x  M2- 

at 


dyW 

at 


The  terms 
have 


^Tsg 

at 


are  generally  negligible  compared  with  the  other  derivatives. 


We  also 


aR(°>  ==  aR(°>  d0g 
3cr  dog  dt 


aR<°) 

de6 


(3.34) 


If  the  station  is  not  moving,  then 


av<°) 

acR 


=0;  otherwise,  this  derivative  depends  on  the 


ground  station  motion  model.  Note  that  although  cR  is  a  receiver  parameter,  it  appears 
in  both  satellite  and  ground  stations  computations  so  that  equation  (3.4)  must  be  used 
rather  than  (3.8)  or  (3.9). 


SUMMARY 

We  summarize  the  parameters  and  the  stages  which  lead  from  them  to  the  doppler 
residuals  d(t)  -  d  t  of  (1.17): 


Satellite  Parameters: 

time,  semimajor  axis,  eccentricity,  inclination,  mean  anomp.iy,  ascending  node, 
argument  of  perigee,  offset  freq. 


<ls  —  (^soj  ao>  eo>  "o,  %  >wo>  fso) 


12 


Ground  Station  Parameters: 
time,  latitude,  longitude,  course,  speed,  frequency 

receiver 

<lG  —  (W  0RO>  ^R0>  ^R0>  SR0>  ^Ro) 
transmitter 

qG  =  (tT0,  $T0,  XT0,  i/T0,  sto,  fo) 


12  Note  that  V  and  U)  as  elements  should  not  be  confused  with  ground  station  course  or  the  earth’s 
angular  velocity. 
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^G^Go) 


Section  4 


qQ(t)  (3-6),  (3.18)  to  (3.20)  R(c)(t^  R(C)(tj  (3^15)  R(S)(tjj  Ris)(t) 

qs(tso)  Appendix  B  ,  r<s)(t),  r<s>(t) 


Hf  (t),  Rg’W 

r<s»(t),  r<s>(t) 

(1.6)  to  (1.9) 

Mfc) 

R(rs)(t),  R$(t) 
r(s)(t),  rCs>(t) 

(1.6)  to  (1.9) 

kuM 

bu<  ^D)  f0,  ^RO> 

(1-14)  to  (1.16) 

.  d(t) 

Derivatives  are  computed  by  the  chain  rule.  The  relevant  equations  are  (3.10)  to  (3.14) 
for  d,  (3.6)  to  (3.9)  for  by  and  bD,  Section  4  for  R'G)  and  R/G)  (note  that  TSG  is  only  a 

linear  factor  in  derivations  with  respect  to  the  above  parameters),  and  Appendix  B  for 
r(S)j  j.(S) 
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4.  MOTION  ON  AN  OBLATE  EARTH 


In  this  section  we  develop  a  four-parameter  model  (i.e.,  initial  latitude,  longitude, 
course,  and  speed)  for  ground  station  motion.  The  obvious  choice  of  great  circle  motion 
is  somewhat  complicated  by  the  consideration  of  a  nonspherical  earth.  Some  of  the  rela¬ 
tions  reduce  considerably  with  the  appropriate  algebraic  manipulation  (e.g.,  (4.18)  and 
(4.19)),  but  the  total  number  of  equations  remains  large.  The  less  interesting  equations 
are  computed  in  Appendix  A  and  summarized  at  the  end  of  this  section.  We  remark 
that,  for  those  only  interested  in  motion  on  an  oblate  earth,  Section  4  and  Appendix  A 
may  be  read  independently  of  the  rest  of  the  report. 

STATION  COORDINATES  TO  GEOCENTRIC  COORDINATES 

The  term  station  coordinates  refers  to  geodetic  latitude,  designated  by  <f>,  longitude 
east,  X,  and  altitude  above  mean  sea  level,  H  (cf.,  (l)  or  Section  3).  Of  course,  since  the 
earth  is  moving,  these  must  also  be  accompanied  by  a  time.  The  other  system,  geocen¬ 
tric  coordinates,  will  be  denoted  by  UK  (  K  =  axis  of  rotation,  I  =  direction  of  vernal 
equinox,  J  =  K  x  I  )  These,  of  course,  also  depend  on  a  specified  epoch.  If  we  let  0g  be 
the  angle  between  the  Greenwich  meridian  and  the  axis  I,  i.e.,  greenwich  sidereal  time, 
then  the  local  sidereal  time  6  is  given  by 

6  =  0t  -1-  X  .  (4.1) 

Let  us  designate  the  equatorial  radius  of  the  earth  by  ae  and  its  eccentricity  by  e,  and  let 


\/l  -  e2sin2<£ 
and 

A  d  q  _  a(e2sin^ 
qs  —  d(sin^)  (l  -  e2sinV)3/2 

Then,  assuming  an  oblate  earth,  the  position  vector  in  IJK  coordinates  is  given  by  [1]  (or 
by  equations  (3.18)  -  (3.20)) 

r(IJK)  __  cos^  y  sjn0,  zy 

where 

x  =  (q+H)  cos 4>  (4.3b) 

z  =  (q  (1-e2)  +  H)  sin<f> ,  (4.3c) 


(4.3a) 


(4.2a) 

(4.2b) 
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and  f  indicates  the  transpose.  Equations  (4.1)  to  (4.3)  describe  the  transformation  from 
station  coordinates  to  geocentric  coordinates. 

Computation  of  the  doppler  shift,  of  the  sensitivity,  or  an  analytic  steepest  descent 
to  estimate  the  state,  requires  the  derivatives  of  these  transformations  with  respect  to 
various  parameters  involved  in  the  station  location  and/or  its  motion.  Since  these 
parameters  are  more  conveniently  expressed  in  terms  of  functions  of  the  station  coordi¬ 
nates  (e.g.,  in  terms  of  sin^  and  cos <f>  rather  than  <f>  itself),  we  retain  that  dependence  in 
computing  the  derivatives.  Let  w  be  an  arbitrary  parameter,  then 


where 


and 


d  Rijk  ,  d  x  „  ,  d  cos 0 
dw  v  dw  dw 


d  x 
dw 


=  qs  cos^ 


d  x  .  a  ,  d  sinfl 

—sm6  +  x— - 

dw  dw 


1±)\ 

dw ' 


d  sin^ 
dw 


+  (q  +  H) 


d  cos  <t> 
dw 


(4.4) 

(4.5a) 


e2)  sin^ 


d  sin# 
dw 


4-  (q(l  -  e2)  4-  H) 


d  sin# 
dw 


q(!  -  e2)  ,  it  d  sin^ 

1  -  e2sin2#  dw 


(4.5b) 


These,  of  course,  do  not  include  the  case  where  w  is  a  function  of  H. 


STATION  MOTION 

We  model  our  motion  as  ’’great  circle”  motion  on  a  nearly  spherical  earth.  By  this 
we  mean  that  the  geodetic  latitude  and  longitude  trace  out  a  great  circle  on  the  unit 
sphere.  Thus,  for  example,  in  Figure  4.1,  P0  is  the  initial  position,  and  P  is  the  position 
after  traveling  the  distance  f  (radians)  along  a  great  circle  specified  by  the  initial  course 
i/0.  This  motion  composed  with  equations  (4.1)  -  (4.3)  is  the  assumed  motion  of  the  sta¬ 
tion.  In  this  sense,  the  local  position  R(t)  is  a  function  of  the  initial  position  tf>0,  Xq>  °f 
the  initial  course  and  of  the  distance  traveled  $(t).  Of  course,  $(t)  is  related  to  the 
initial  speed  s0  and  the  time  elapsed  t  -  t0  as  diagramed  below. 


to,  K  Xo,  J'o,  So  to,  K  »0,  it)  (4'4Q)"(4,4|)  mm  (4-lH4-3),  R(t) 


For  lack  of  a  better  term,  we  call  this  ’’pseudo  great  circle  motion.” 

The  velocity  R  is  then,  of  course,  the  time  derivative  of  this  motion.  It  will  sim¬ 
plify  our  results  to  separate  the  part  of  the  motion  due  to  the  earth’s  rotation  from  that 
of  the  station  relative  to  the  earth.  Let  Cl  be  a  vector  along  the  earth’s  axis  with  magni¬ 
tude  equal  to  the  earth’s  rate  of  rotation  (cf.,  (3.5))  and  v  the  velocity  the  station  would 
have  if  we  instantaneously  stopped  the  earth’s  rotation.  Then 
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R(t)  =  fi(t)  x  R(t)  4-  v(t) . 


(4.6) 


The  vector  fi(t)  x  R(t)  is  simply 


dR  dfl  13 

de  dt‘ 


Hence,  v  is  computed  by  taking  the  time 


derivative  of  R  but  setting  d#/dt  =  0. 

One  must,  however,  resist  the  temptation  to  further  simplify  by  trying  to  directly 
convert  the  course  y  and  speed  s  to  rectangular  coordinates  (e.g.,  north  and  east)  and 
rotating  that  to  IJK  coordinates,  thus,  apparently  obtaining  v  without  the  need  to 
differentiate  (4.2)  to  (4.3).  More  precisely,  let  the  transformation  from  SEZ  (topocentric: 
south,  east,  up  )  coordinates  to  IJK  be  given  by  the  rotation  ((!]  or  equation  (3.22)) 


where 


V(UK)  =  d-iv(sez) 


D-1  = 


sin0  cos#  -sin# 
sin^  sin#  cos# 
-cos  ^  0 


COS<£  cos # 
cos <f>  sin# 
sin^ 


(4.7) 


(4.8) 


and  is  valid  for  an  ellipsoidal  earth.  One  might  expect  the  velocity  vector  to  lie  in  the 
SE  plane  at  an  angle  of  ir  -  v  to  the  south,  i.e., 


and 


V(SEZ)  _£.  £  aR 


-cosy 

siny 

0 


V(IJK) 


j  aR  D"1 


-cosy 
siny  , 
0 


(4.9) 


(4.10) 


where  j  indicates  the  time  derivative  of  and  aR  is  the  radius  of  the  earth  at  R. 

However,  this  is  only  approximately  correct.  Due  to  the  earth’s  oblateness,  the 
actual  curve  traced  on  the  earth’s  surface  is  not  the  great  circle  of  a  sphere,  and  the  SEZ 
velocity  vector  is  only  approximately  given  by  (4.9).  The  exact  equations  are  computed 
below.  We  also  note  that  we  must  specify  whether  $  is  a  constant  or  whether  the  speed, 
||  v  ||  is  a  constant.  Although  these  are  approximately  equivalent,  once  again,  because 
of  the  oblateness  of  the  earth,  they  are  not  exactly  the  same.  The  final  result  is 


v(;JK)  =  j  D-1 


a 

-cosy 

siny 

+  b 

cos <t>  cosy 

0 

' 

0 

0 

j 

(4.11) 


13  fi  depends  on  t  since  the  earth’s  axis  precesses.  fi,  which  is  negligible,  is  effectively  set  to  zero  if 
we  use  equation  (4,6);  however,  if  one  is  fussy,  one  may  replace  the  first  term  in  (4.6)  with 

dR  d  0 
do  dt 
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■  xr*  xfws  r/v  x\i  xrs/  v/vvtv  *sr\  v* 


ivitvi  Armn  atma-a/i  rvrt  mt  *xi  it/vcrt  n  n  Kr.-vn-ijt  vji  va  vkv 


where 


and 


a  =  q  4-  H, 


(4.12a) 


■^-4-  cos^  =  <ls 

V 


COS  <j> 

sin^ 


(4.12b) 


We  proceed  to  derive  (4.11).  Using  definitions  (4.12),  equations  (4.3)  and  (4.5)  may 
be  rewritten 


x  =  a  cos <j> 

(4.13a) 

z  =  a  s\n<f>  -  e2  q  sin^ 

(4.13b) 

dx  d  cos<j>  .  ,  .  ,  d  sin<£ 

=  a  +  b  sin^ 

3w  ow  a  w 

(4.13c) 

9  z  =  a  9  sin^  b  cosd,  d  sin^ 

3w  aw  9  dw 

(4.13d) 

In  order  to  determine  the  velocity,  we  wish  to  set  w  =  £  and  use  v 
Equations  (A. 15a)  and  (A. 15b)  give 

=  (3R/3f)  f  • 

d  sin^  , 

- - — —  =  COS 0  cost' 

a? 

(4.14a) 

d  cosd>  •  , 

— — — -  =  -  sin  <p  cost /, 

(4.14b) 

while  from  (A.38a)  and  (A.38b) 

3sin0  =cos6  sint' 
di ;  cos  <j> 

(4.15a) 

and 

5  cos 9  ■  a  sinf 

— - —  =  -  smv - - 

a?  cos  <p 

(4.15b) 
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Combining  (4.13c)  and  (4.13d)  with  (4.14a)  and  (4.14b)  yields 


^  X  ■  • 

— —  =  -a  sin^  cosy  -f  b  sin$  cos <f>  costs 
d$ 


(4.16a) 


=  a  cos  <f>  cos  is  -  b  cosV  cosy 

3? 


(4.16b) 


Substituting  (4.15)  -  (4.16)  into  (4.4),  we  obtain 


d  r(1j1<) 

a? 


-a  sin$  cosy  cose?  -  a  sine?  siny 

b  sin^  cos  cj>  cos  is  cos  <? 

-a  sin^  cosy  sin(?  +  v.  cos 0  siny 

+ 

b  sin^  cosS  cosy  sine? 

a  cos <j>  cosy 

-b  cos 2^  cos  y 

(4.17) 


Comparing  this  with  (4.8),  we  find 

a  r(ijk> 

a? 

Since 


-cosy 

cos <}>  cosy 

a 

siny 

+  b 

0 

0 

0 

d 

v  =  f 


a  r 

a?  ’ 


(4.18) 


(4.19) 


this  proves  (4.11). 

It  remains  to  describe  the  relationship  of  £  to  speed.  Since  speed  is  the  magnitude 
of  the  velocity,  we  have 


where 


i  =  speed/  ||  || 


(4.20) 


||  || 2  =  a2  +  b2  cos '<{>  cos2y  -  2  a  b  cos  <£  cos' is 

d<; 


(4.21) 


from  (4.18)  and  the  orthogonality  of  the  matrix  D.  Thus,  we  see  that  we  may  postulate 
either  motion  with  a  constant  speed,  or  with  a  constant  rate  along  the  great  circle  (i.e., 
constant  j),  but  not  both.  Although  speed  is  almost  always  the  measured  parameter, 
constant-speed  motion  would  require  a  complex  integration  of  the  function 
|j  (dR/Sj)-1  ||  along  the  great  circle  path  in  order  to  determine  f  as  a  function  of  time 
(the  integral  of  (4.20)  ).  Consequently,  it  seems  advisable  to  model  the  motion  by 


?  =  s0  (t  -  t0) 


(4.22) 
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where  s0  is  a  constant  set  equal  to  the  true  j  at  time  t0.  In  general,  this  is  an  excellent 
approximation  (it  would  be  exact  for  a  spherical  earth),  and  the  error  involved  is  less 
than  that  involved  in  the  assumption  of  pseudo  great  circle  motion.  (In  most  cases,  the 
actual  motion  is  only  approximated  by  a  set  of  constant-course  segments.  Even  for  great 
circle  motion  on  a  spherical  earth,  the  course  is  always  changing.)  Furthermore,  for 
computation  of  the  doppler,  we  may  still  use  the  exact  value  of  j.  More  precisely,  at  a 
constant  speed  spd,  (4.20)  at  time  t0  yields 


S.  —  j(t„)  =  —gp -  , 

ii  f^to)  ii 

(4.23) 

and  at  time  t,  we  have 

.  II  #(*o>  II 

m\—  sPd  _  -  of 

^  ati  u  sb 

II  ^(t)  II  II  ^(t)  || 

(4.24) 

In  summary,  we  propose  two  possible  models,  both  of  which  avoid  the  integration 
of  (4.20)  with  respect  to  time.  In  the  first  model,  we  assume  that  f  =  s0,  a  constant. 

f  (t)  =  s0(t  -  t0) 

fiV 

(4.25a) 

UJ- 

j(t)  =  s0. 

(4.25b) 

This  model  is  physically  consistent,  i.e.,  j  is  the  time  derivative  of  f.  However,  it  is  clear 
from  (4.20)  that  in  this  case  the  speed  will  not  be  a  constant.  For  the  second  model,  we 
define  f  and  £  by 


(ii): 


f  (t)  =  s0(t  -  t0) 


11  ffito)  11 

II  frMII 


(4.26a) 


(4.26b) 


In  this  case  £  is  not  the  time  derivative  of  £,  and  v  of  (4.19)  will  not  be  the  velocity 
corresponding  to  the  motion  (4.26a).  However,  v  will  be  the  velocity  at  the  position  $(t) 

3D 

given  that  the  speed  is  s0  ||  (t0)  || .  Thus,  if  the  true  motion  is  constant  speed, 

of 

model  (ii)  will  give  a  (slightly)  incorrect  position  but  a  correct  velocity  v  for  that  posi¬ 
tion. 
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DERIVATIVES  OF  STATION  MOTION 

From  (4.19)  we  have,  for  an  arbitrary  parameter  w, 

dv  •  d  dR  ,  d$  dR 

3w  ^  3w  d$  dw  d$ 

If  w  =  s0,  then  in  both  models  (i)  and  (ii)  equation  (4.27)  becomes 

dv  •  d  dR  j_  dR 

ds0  ?  ds0  d?  s0  df 


(4.27) 


■n  >\  52R  . 
=  ?(t-t0)—  + 


so 


dR 

a? 


(4.28) 


For  model  (i)  these  simplify  to 

(i):  J^=s»ffr  for'V74s°  (4'29a) 

dv  /.  .  %  a2R  ,  dR  9Q.  \ 

(,):  -^=s°(t-to)-ar+ir  (  ’ 


On  the  other  hand,  model  (ii)  requires  considerably  more  computation.  Since,  from 
(4.24), 


iL  =  _ 

aw 


iif  ii 


a 

aw 


aR 

3? 


we  have  for  model  (ii) 

(ii): 


dv 

aw 


=  $ 


a  aR 


a w  a? 


aR 

a? 


aR 

3? 


-i  A_ 

aw 


(4.30) 


Ilf- II 


(4.31a) 


(ii): 


dv 

dSn 


=  ? 


fi  t  \  32R  .  1  dR 


(4.31b) 


To  a  very  high  approximation,  we  could  drop  the  second  term  in  (4.31a);  however,  if  we 
use  the  partials  in  a  steepest  descent,  then  convergence  to  the  true  minimum  of  a  cost 
function  can  only  be  guaranteed  if  the  second  term  is  retained. 
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Let  us  derive  a  slightly  more  explicit  form  for  — —  (dR/d$).  From  (4.18), 


d  dR 
dw  d$ 


(-a  4-  b  cos^) 
a  siny 
0 


cosy 


(4.32) 


Definitions  (4.12)  and  (4.2)  give  us 

e2  o  o 

-a  +  b  cos^  =  -  q  -  H  4 - -  q  cos V 

ae” 


=  -  II  4-  q  (-14- 


1  -  e2  sin  i<j) 


cos  2<j>) 


=  -  H  - 


(4.33) 


Then,  indicating  3  /3w  by  subscripting  w,  we  have 


d_ 

c)w 


(-a  4-  b  cos^)  cosy 
a  siny 
0 


4-  D-1 


(-a  4-  b  cos^)w  cosy  4-  (-a  4-  b  cos^)  (cosy)w 
aw  siny  4-  a  (siny)w 
0 


(4.34a) 


where,  from  (4.12a)  and  (4.33), 


aW  -  9s 


3sin<ft 

3w 


(4.34b) 


(-a  4-  b  cos0)w  = 


3(1  -  e2)  2  dsimft 

,2  4  4s  aw  ' 


(4.34c) 


SUITABLE  STATE  VARIABLES  FOR  STEEPEST  DESCENT 

In  many  cases,  the  gradients  computed  in  the  previous  section  will  be  used  in  a 
steepest  descent  or  sensitivity  calculation  for  which  the  state  variables  include  course 
and  speed.  It  can  be  be  seen  from  (4.19),  (4.22)  and  the  fact  that  sinf  appears  as  a 
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factor  in  dR/di/0  (substitute  successively  (4.57),  (4.50),  (4.46),  and  (4.5)  into  (4.4)),  that 
at  Sg  =  0,  the  partials  of  R  and  v  with  respect  to  are  zero,  and  the  resulting  ”A” 
matrix  (cf.,  (2.10))  will  be  singular.  This  is  essentially  an  artifact  of  the  ’’polar”  (i.e., 
course,  speed)  coordinate  system  and  may  be  alleviated  by  changing  state  variables. 

A  natural  choice  for  new  state  variables  is  the  rectangular  velocity  coordinates 
defined  by 


Pno  =  s0  COSI/q 


(4.35) 


Pt o  ==  s0  sini/Q 


1 


where  the  subscripts  n  and  e  are  intended  to  indicate  velocity  north  and  east  respec¬ 
tively.  We  can  transform  to  this  system  using  the  Jacobian 

d&o  dPto 
ds0  di/0 

or,  more  practically,  its  inverse 


COSl^o 
Sg  sin  i/0 


smi/0 

Sq  COSI/q 


(4.36) 


ds0  dvQ 

dpM  dp  to 


cosi/0 

sini/0 


sini/0 

s0 

cosi/0 

s0 


(4.37) 


For  example,  from  the  chain  rule,  (4.29)  and  (4.37)  yield 


0): 


’  av 

sini/Q 

COSI/q 

— - 

dPno 

So 

dv 

COSI/Q 

dPa  . 

sint^o 

So 

,  v  a2R 

so(t  - 


+ 


OR 


8  dR 
dv0  a? 


cosi/q  -  sint'o 
sinz/0  cos^q 


So(t  ~  *0) 


a2R 

a? 


+ 


aR 

5? 


a  aR 

a^0  d$ 


cosu0  -  sini^o 
sinj/0  cost'o 


a2R  ,  aR 

a  9R 

di/0  a? 


(4.38) 
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The  same  coordinate  changes  circumvent  the  vanishing  of  dR/du0  at  s0  =  0.  To 
see  this,  we  first  note  that  the  partials  of  R  with  respect  to  /?no  and  /?eo  take  the  form  of 
(4.38)  with  the  vector  on  the  right-hand  side  having  components  dR/ds0  = 
(t  -  t0)  dR/d?  and  (l/s0)dR/dy0  ,  i.e.. 


'  <9R 

r(t  -  to)  iR  1 

3A>0 

COSl^g  -  SUli/g 

v  01  a? 

dR 

sint/g  COSI/g 

JL 

.dfito  . 

So  dl/g 

(4.39) 


From  (4.4)  and  (4.5),  we  observe  that  every  term  of  dR/di>0  contains  one  of  the  four 
expressions  ds\n<j>/duQ,  dcos^>/du0l  dsin 0/duo,  or  dcos8/dis0.  From  (4.46),  (4.50),  and 
(4.57),  we  also  see  that  each  of  these  contains  a  factor  of  sin?.  Thus,  (l/s0)  dR/du0  con¬ 
tains  a  factor  of  sin?/sg  which  has  the  finite  (and  generally  nonzero)  limit  of  (t  -  to)  as  s0 
goes  to  zero. 

Another  aspect  which  should  be  mentioned  is  the  potential  for  negative  ?.  It  is 
entirely  possible,  during  a  steepest  descent  procedure,  to  drive  the  state  variable  so  nega¬ 
tive.  This  poses  no  problems  provided  we  interpret  it  as  motion  in  the  opposite  direc¬ 
tion  to  the  indicated  course.  More  precisely,  suppose  we  are  faced  with  a  negative  ?, 
and,  hence,  negative  ?.  We  may  define  primed  variables  by  ?'=-?; 
i/0  =  i/0  +  7r;  and  iJ  —  v  +  ir.  Equations  (4.40)  to  (4.43)  of  the  next  subsection  remain 
invariant  under  this  change  of  variables.  Since  f  is  positive,  we  may  interpret  the  nega¬ 
tive  ?  to  be  the  same  motion  as  -?  with  the  course  reversed,  i.e.,  with  the  course  replaced 
by  (course  4-  7r).  More  simply  stated,  we  may  work  with  equations  (4.40)  to  (4.58)  with 
impunity  even  if  ?  is  negative.  Also,  if  a  particular  set  of  computations  results  in  nega¬ 
tive  ?,  then  the  course  u0  should  be  interpreted  as  the  direction  opposite  to  the  physical 
course. 


SUMMARY  OF  SPHERICAL  TRIANGLE  EQUATIONS 

For  convenience,  we  assemble  here  the  equations  from  Appendix  A,  which  are 
needed  to  compute  sin^,  cos<f>,  sinu,  cost',  sin#,  cos#,  and  their  derivatives.  Note  that 
when  unsubscripted  variables  appear  on  the  right-hand  side,  they  must  be  computed  via 
previous  equations. 


sin^  =  sin^g  cos?  +  cos^g  sin?  cosfg 

(4.40a) 

cos^  =  \/l  -  sin  V  <t>  e  (— ?r/2,  ir/2) 

(4.40b) 

cos?  -  sin^g  sin^ 

COSA  =  - - - 7 — 

COS 00  COS0 

(4.41a) 

sinfsinz/0 

sinA  = - - — 

cos  <p 

(4.41b) 

costs  =  cost'g  cosa  -  sint'g  sinA  sin^g 

(4.42a) 

sim/g 

sin^  =  cosipo - — 

cos  <p 

(4.42b) 
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sin0  =  cos(X0  +  0g)  sinA  +  sin(X0  +  0g)  cosa 
cos 6  =  cos(X0  4-  0g)  cosa  -  sin(X0  +  0g)  sinA 


(4.43a) 

(4.43b) 


=  cos?  cos#0  -  sin?  cosi/q  sin#0 
d#  o 

dcos #  _  sin#  dsin # 

d<j>0  cos  #  d#  0 

dsin#  _  dcos #  _ 

d\0  ~  dx0 

dsin#  ,  .  . 

--■■■  =  -cos#0  sin?  siny0 

OUq 

dcos #  _  sin#  dsin# 

dy0  cos#  dy0 

=  cos#  cosy 

d? 

dcos  #  .  , 

— - — —  —  -sin©  cosy 

a? 


(4.44a) 

(4.44b) 

(4.45) 

(4.46a) 

(4.46b) 

(4.47a) 

(4.47b) 


dsinA  _  sinA  dcos# 
d#  o  cos#  d#0 

dcosA  __  sin#  +  cosa  gjn^  _  sinffo  dsin# 
d# o  cos#  cos#0  cos#0  cos#  d#0 

cosa  dcos# 
cos#  d#0 

dsinA  _  dcosA  _ q 

dXo  dXo 

dsinA  __  sin?  cosy0  sinA  dcos# 
dy0  cos#  cos#  dy0 

dcosA  _  sin#0  dsin#  _  cosa  dcos# 

dy0  cos#0  cos#  dy0  cos#  dy0 


(4.48a) 

(4.48b) 

(4.49) 

(4.50a) 

(4.50b) 


dcosy 

d#0 


dsiny  _  sin#0  siny0  siny  dcos# 

d#  o  cos#  cos#  d#0 

dcosA  •  dsinA  .  .  .  , 

c°sy0  — - siny0  — — —  sin#0  -  smy0  smA  cos#0 

d#0  d#  o 

dsiny  _  dcosy  __  q 
dX0  dX0 

dsiny  _  ,  cost/o  siny  dcos# 

dy0  cos#  cos#  dy0 


(4.51a) 

(4.51b) 

(4.52) 

(4.53a) 


58 


dcosy 

du0 


.  ,  3cOSA  .  .  , 

-sin  i^o  cosa  +  cos z^o  — - cos^o  sin  A  sin<po 


dun 


Ssiny 

a? 

dcosy 

a? 


-  siny0  sin? 
sin<f> 


dsinA 
dv. 


o 


i  siny  cosy 


COS0 


.  2  sin<i 
=  -  snry  — 


cos<f> 


(4.53b) 

(4.54a) 

(4.54b) 


dsin#  A  ,  0  <.  dsinA  ,  .  a  \  <9cosa 

_  _  COS{X0  +  *.)  -5£-  +  ™(>.  +  »,) 

Ozsl  _  cos(X0  +  <y  -  sin(X0  +  «g) 

dsinfl 

ax0 

dcos0 

axn 


dfa  d<i>0 

=  -sin(X0  4-  9g)  sinA  4-  cos(X0  +  9g)  cosa 

=  -sin(X0  +  9g)  cosa  -  cos(X0  +  9g)  sinA 


3sin0 

5y0 

=  cos(X0  +  0g) 

dsinA 

ay0 

+  sin(X0 

9cos0 

ay0 

—  cos(X0  +  0g) 

dcosA 

ay0 

-  sin(X0 

dsinO 

=  cos  9 

siny 

a? 

cos  cj> 

dcos  9 

=  -  sin*  Sini/ 

a? 

cos  <p 

dcOSA 

ay0 

dsinA 

ay0 


(4.55a) 

(4.55b) 

(4.56a) 

(4.56b) 

(4.57a) 

(4.57b) 

(4.58a) 

(4.58b) 
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Appendix  A:  GREAT  CIRCLE  MOTION  AND  ITS  DERIVATIVES 


A  standard  spherical  triangle  is  pictured  in  Figure  A.l.  We  are  considering  motion 
along  the  great  circle  arc  £  from  Pj  to  P2  where  the  latitude  fa  and  course  i>\  of  the  ini¬ 
tial  point  as  well  as  the  distance  f  (in  radians)  an  known.  (The  longitude  Xj  is  also 
known,  but  not  relevant  at  the  moment.)  We  wish  ti>  determine  the  new  latitude  fa,  the 
change  in  longitude  A,  the  new  course  i/2)  and  their  partial  derivatives.  Figure  A.l,  as 
well  as  standard  spherical  trigonometric  formulae,  assumes  that  all  angles  and  sides  of 
the  triangle  in  question  lie  in  the  interval  [0,  tt).  This  is  true  of  the  colatitudes  7t/2  -  fa 
and  7r/2  -  fa ;  however,  uh  f,  a,  and  v2  will  generally  be  in  the  interval  0  to  2n.  To  han¬ 
dle  this  problem,  we  transform  to  variables  which  satisfy  the  desired  conditions  and  dis¬ 
tinguish  them  by  the  addition  of  a  prime.  Four  cases  occur,  as  illustrated  in  Figure  A.2 
and  summarized  in  Table  A.l. 


Table  A.l.  Transformations  to  standard  spherical  triangle. 


For  those  transformations  omitted,  the  primed  value  equals  the  unprimed 


value. 


We  could  perform  all  our  derivations  assuming  case  (0)  and  then  justify  the  results  for 
all  cases  by  an  appeal  to  analyticity  of  the  individual  variables  (which  implies  a  unique 
extension).  However,  with  a  bit  more  labor,  we  achieve  the  same  end  directly,  using  the 
following  sign  relationships  (easily  verified  from  Table  A.l): 

cosa  =  cosa'  (A-l) 

cosf  =  cos^  (A.2) 


sin? 

sin^ 


cosiq  cos/A) 

cos  cosiA 


(A.3) 


A-l 


sini/j  sin  1^2 

sin  sinz/2 

sinA  =  sin?  sin^ 
sinA'  sin?'  sini/,1 

Equations  (A.3)  and  (A.5)  imply 

cosi/2  sini/j  sinA 

cost/1 o  sini^!  sinA' 

Table  A.2  contains  a  summary  of  the  locations  of  the  various  equations. 


(A.4) 


(A.5) 


(A.6) 


Table  A.2.  Equation  numbers  for  various  computations  . 


indep 

variable 

sin/cos 

d/dfa 

erivatives 

a/ax 

of  sin/cos 
d/dvx 

a/a? 

<t>2 

A.7 

A.10 

A.ll 

A. 13 

A. 15 

A 

A.8 

A.20 

A.21 

A.22 

Vo 

A.9 

A.25 

A.26 

A.27 

A.30 

Oo 

A.32 

A.33 

A.34 

A.35 

A.  38 

FUNCTIONS  OF  P2 

In  the  equations  which  follow,  we  shall  assume  that  no  point  of  the  arc  passes 
through  one  of  the  poles.  To  solve  for  >  we  use  the  law  of  cosines'  applied  to  the  sides 
of  the  spherical  triangles  in  Figure  A.2.  Thus, 

sin^2  =  sin^!  cosj7  +  cos^j  sin^  cosi^ 


=  sin^j  cos?  +  cos^i  sin?  costq, 


(A.7a) 


and 


f  Our  choice  of  equations  in  what  follows  is  dictated  by  a  desire  to  avoid  indeterminate  forms  and 
simplify  sign  determination.  Thus,  only  cosines  of  the  latitude  appear  in  denominators. 
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(A.7b) 


COS  1^2  =  %/ 1  —  sin2^2  <f> 2  €  (— 7r/2,  7t/2). 

Note  that  we  have  used  (A.2)  and  (A.3)  in  the  derivation  of  (A.7a). 
Similarly  for  a', 


cosa'  = 


which,  by  (A.l)  and  (A.2),  implies 


cosa  = 


cos^-  sin^j  sin02 
COS^j  COS^2 


cosf  -  sin^!  sin^2 

COS^i  COS (^2 


From  the  law  of  sines, 

sin^sini/j 

sin  a'  = - - — 

COS02 


so  that  (c.f.,  (A.5)) 

sinf  siniq 

sinA  - - - — 

cos  92 


(A.8a) 


(A.8b) 


The  law  of  cosines  for  angles  yields 

cosi/2  =  -cos cosa1  +  sinj/j  sinA*  sin^; 


that  is,  by  (A.l),  (A.3),  and  (A.6), 

cos i^2  =  cosi^j  cosa  -  sini/j  sinA  sin^j. 
Again,  from  the  law  of  sines, 


,  ^  sini/i 

smi/2  =  cos0, - — , 

cosp2 

or,  using  (A.4), 

sini/j 

sini'2  =  cosd>i - — , 

cos^2 


(A.9a) 


(A.9b) 


Note  that  in  this  section,  as  in  those  which  follow,  some  expressions  may  contain  quanti¬ 
ties  which  must  be  computed  from  previous  equations.  For  example,  (A.9b)  depends  on 
cos^2,  which  may  be  computed  using  (A.7). 
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DERIVATIVES  OF  sin^o  and  cos^2 
From  equations  (A.7a), 


3sin^2 

— - - =  cosf  cos<f>i  -  sinf  cosi/j  sin^j 

o<Pi 


(A.  10a) 


Also, 


where  we  have  used 


dcos<j>2 

sin$2  9sin^2 

(A.  10b) 

dh 

C0S^2  9<j)  i 

3sin<£2 

3cos^2 

~0  ’ 

(A.U) 

3cOS^2 

a... 

s\n<f>2  9sin^2 

_ j.  a...  > 

(A.12) 

3w  cos^2  3w  ’ 


which  follows  from  (A.7b)  for  arbitrary  w. 

The  partials  with  respect  to  Vy  are,  from  (A.7a)  and  (A.12), 


dsm^2 

— - - =  -cos^i  smj  sin^i 


(A.  13  a) 


dcosfio  sin ^2  9sin^2 

diq  cos^2 


(A.  13b) 


Before  computing  the  derivatives  with  respect  to  £,  we  note,  from  Figure  A.l,  that 
taking  the  derivative  with  respect  to  ?  at  P2  is  the  same  as  taking  the  derivative  at 
P2  =  Pj  by  setting  $■  equal  to  zero,  and  then  replacing  the  subscript  2  with  the  subscript 
1.  That  is,  the  rate  of  change  of  sin^  or  cos^  with  respect  to  £  at  an  arbitrary  point 
depends  only  on  that  point.  As  we  shall  see,  this  provides  us  with  a  quick  method  for 
simplifying  the  resulting  expressions.  Following  the  indicated  procedure,  we  have  from 
(A.7a) 


3sin^2 

d<; 


(-sin^i  sin£  +  cos^  cosj  cos^j) 


A- 6 


=  COS^i  COSiq. 


(A.14) 


Then,  since  (A.14)  is  only  a  function  of  Pt  and  since  Pi  in  this  context  is  arbitrary,  we 
may  replace  it  by  P2  to  get 


Using  (A.12),  we  find 


dsin<^2 

a? 


=  COS02  COSJ/2. 


dcos^o 

— — —  =  -sin^2  cosu2  . 


(A.  15a) 


(A.15b) 


For  completeness,  let  us  show  that  this  is  equivalent  to  just  taking  the  derivative  of 
(A.7)  with  respect  to  We  have  from  (A.8a) 


cos^2  = 


cosj  -  sin^i  sin^2 


so  that 


COSI/j  cosa  cos^2  =  COSt/j 


COSA  cos^i 

cosf-  sin^i  sin^2 

COS^i 


(A.16) 


- (cosf  -  sin2^j  cosf-  sin^j  cos^j  sin£  cosi/J 


COS^!  COS?  2  2  ■  1  ■ 

= - 7 - cos  <j>i  -  cos  swupi  smf  ,  (A.17) 

COS  <Pi 

where  the  second  last  line  follows  from  substitution  of  (A.7a).  Also,  from  the  law  of 
sines  for  Figure  A.l  (equation  (A.8b)), 

cos^2  sinA  sini/j  sin^  =  sinj  sin2^  sin^.  (A.18) 

Multiplying  (A.9a)  by  cos^2  and  substituting  (A. 17)  and  (A.18)  in  the  resulting  expres¬ 
sion,  we  get 


cos ^2  cosi/2  =  cos cos^j  cosj-sin^j  sinf 


(A.19) 
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Then  (A.15a)  equals  the  left-hand  side  of  (A. 19)  while  the  right-hand  side  is,  indeed,  the 
derivative  of  (A.7a)  with  respect  to 

DERIVATIVES  OF  sinA  AND  cosa 
Differentiating  (A.8b),  we  have 


dsinA  sinA  dcos^ 

d(j>i  COS^j  9<f>  i 


From  (A.8a) 


3cosa 

d<f>i 


sin 

COS<f>2 


+ 


COSA 

COS^! 


sin^!  - 


sin^x 

COS^i  COS^2 


Ssin^2 

~w 


cosa  doosif)^ 
COS  1^2  d<f>i 


Finally, 


3sinA  _  5cosa 

aXj  “  d\t 


(A.20a) 


(A.20b) 


(A.21) 


The  partials  v.  ith  respect  to  iq  are  also  found  by  differentiating  equations  (A.8). 
From  (A.8b) 


5sinA  _  s*n?  cosi^i  sinA  dcos<j> i 

dtq  cos^2  cos(j>o  di/i 

Similarly,  from  (A.8a), 

dcosA  sin^j  5sin^o  cosa  dcos<f>n . 

di>i  cos^i  cos^2  cos^2  diq  ’ 


(A.22a) 


(A.22b) 


For  $•  we  find  that  we  shall  only  need  the  derivatives  at  f  =  0,  so  we  proceed  as  in 
the  derivation  of  (A. 14).  Equations  (A. 8a)  and  (A. 15)  imply 


5cosa 

J  sin? 

COSA 

dcOS(j>2 

sin^q 

asin^2  1 

a? 

0 

1  COS^i  COS^2 

COS^2 

* 

COS^j  COS<f>2 

a?  J 

=  0  - 


COS  <f>i 


(-sin^j  costq)  - 


(cos^i  cosiq) 

COS' <P\ 
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=  0. 


(A.23) 


From  (A.8b) 


dsinA 

3? 


cos?  sini/j 
COS  ^2 


sinA 

COS<fi2 


3cOS^2  ^ 
~dl~' 


0 


sini/j 

COS^i 


(A.24) 


DERIVATIVES  OF  sinj/2  AND  cosi^ 
From  (A.9b) 


dsini^  sin^j  sini/j  sint/2  dcos <f>2 
d<f>i  cos^2  costf>2  d<f>i  1 


(A.25a) 


and  from  (A.9a) 


dcosu2 

~dh~ 


COSI^j 


dcosA 

d<t>i 


■  -  sini/i 


~^nA  sin^j  -  sini^  sinA  cos <f>x  . 
d<P  i 


(A.25b) 


Of  course, 


3sinv2 

d\x 


dcos u2 
3Xj 


(A.26) 


Next,  taking  derivatives  of  (A. 9)  with  respect  to  ux,  we  get 


9sin^2 

— - - =  COS?)! 

OUi 


COSl'i 

cos<f>2 


sin  v2  dcos<f>2 
COS  ^2  d^l  ’ 


and 


dcosu2  .  0COSA  .  . 

- =  -siniA  cosa  +  cosi/j  — - - cosi/j  sinA  smpi 

aux  ovx 


-  sini^  sin^i 


dsinA 

dvx 


(A.27a) 


(A.27b) 
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We  employ  the  usual  trick  to  obtain  the  partials  with  respect  to  g. 
and  (A.15b): 


From  (A.9b) 


3sini/o 

dg 


sin  1/2  3cos^2 
cos<j>2  dg 


0 


siniq 

- —  simii  cost/,. 

COSipy 


Similarly,  from  (A.9a)  and  (A.24), 


Scos/'o 

dg 


Scosa  .  .  ,  SsinA 

cos Vy  — — - sini/j  sirupi  — — — 


=  -  simq  sin^j 


siniq 

COS  <j>i 


Thus, 


dsini/o 

— - — -  =  sin  1/0  cosi'o 
dg 


sin^ 

cos^2 


dcosi/2 

~eT 


=  -  sin 2i/2 


sin^2 

COS^2 


(A.28) 


(A.29) 


(A.30a) 


(A.30b) 


DERIVATIVES  OF  sin02  AND  cos02 

Instead  of  functions  of  the  longitude,  we  are  generally  concerned  with  those  of  6, 
the  sidereal  time, 


#2  —  +  Xj  +  A  ,  (A.31) 

where  Xj  +  A  is  the  longitude  at  P2,  and  0g2  is  the  Greenwich  sidereal  time  at  point  P2  , 
i.e.,  when  the  station  is  at  P2.  Equation  (A.31)  implies 


sin02  =  cos(X!  4-  0g2)  sinA  +  sin^  +  0g2)  cosa 


(A.32a) 


Thus, 


cos02  =  cos(Xj  +  0ga)  cosa  -  sin(X!  +  0ga)  sinA  (A.32b) 


dsin62  /N  ,  „  ^  SsinA  ,  a  \  Scosa 

--  =  cosfx,  +  <y  + «»(>■.  +  ej  — 


(A.33a) 


dcos$2  3cosa  .  A  ,  Q  \  SsinA 

_  =  costX,  +  IJ  -gj-  -  s,n(X,  +  »J  -  jjj- 


(A.33b) 


Ssin^-i 


-  =  -sin(X!  +  0ga)  sin  A  +  cos^  +  0ga)  cosa  (A.34a) 


dcos#2 


=  -sin(Xj  +  0ga)  cosa  -  cos(Xj  +  0ga)  sinA  (A.34b) 


dsm02  SsinA  ,  .  ,x  ,  „  \  3cosa 

— —  =  cosfX,  +  <y  -rr—  +  sm(X,  +  «J  -j— 


(A.35a) 


dcosd2  .  .  5cosA  .  /x  ,  /i  \  3sinA 

—  =  cos(X,  +  9J  -jjp  -  -  A.  +  -ST- 


(A.35b) 


From  (A.32),  (A.23),  and  (A.24),  we  have 


dsin02  ,n  ^  SsinA  ,  n 

— -  =  cos(02  -  a)  —t—  +  0 

d$  A=0  <0  A=0 


=  COS#! 


COS^i  ’ 


(A.36) 


and  similarly 


dcos#i  .  „  sint^ 
— - —  =  -  sm0! - — 

COSfj 


(A.37) 


A-ll 


Again,  since  this  must  hold  for  all  P*  along  the  path^,  we  have 


dsinflo 


smi/2 

=  COS»2  - T~ 

COS^2 


(A.38a) 


dcos02 

~di~ 


sini/2 

=  -  sin02 - — . 

COS^2 


(A.38b) 


t  Note  that  the  definition  of  8  as  a  function  of  P  is  independent  of  the  point  P1(  which  permits  re¬ 
placing  the  subscript  1  with  2  in  (A.36)  and  (A.37).  This  was  not  possible  in  (A.23)  or  (A.24)  since 
A  is  defined  in  terms  of  Pj.  Of  course  ft  in  (A.36)  also  depends  on  P1(  but  if  we  define  ft,  as  the  arc 
length  relative  to  another  point  Pb,  then  we  have  ft,  =  ft-  arclength  PiPb  =  ft  -  constant,  so  that 
d/dft  =  d/dftb,  which  is  independent  of  the  choice  of  the  point  Pb. 
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Appendix  B:  KEPLER  MODEL  AND  DERIVATIVES 


The  Kepler  model  is  specified  by  six  elements  and  an  epoch: 


a  =  orbit  semimajor  axis 
e  =  orbit  eccentricity 
M0  =  mean  anomaly  at  epoch 
i  =  inclination 

Cl  =  longitude  of  ascending  node 
oj  —  argument  of  perigee 
t0  =  epoch  of  elements. 


The  first  three  quantities  determine  the  position  and  velocity  for  all  time 

with  respect  to  the  orbital  plane  and  relative  to  perigee,  i.e.,  in  the  so-called  perifocal 
coordinate  system  (P  is  the  semimajor  axis  to  the  periapsis,  Q  rotated  90°  in  the  direc¬ 
tion  of  motion,  and  W  determined  by  the  right-hand  rule  (l|).  The  remaining  three  ele¬ 
ments  determine  the  orientation  of  the  orbital  plane.  We  then  have,  with  IJK  represent¬ 
ing  geocentric  coordinates  at  epoch  t0, 


r(lJK)  =  u  r(PQW) 


j(UK)  =  u  *(PQW) ; 


(B.l) 


where 
U  = 


cos  Cl  cos  oj  —  sin  Cl  sin  oj  cos  i  -cos  Cl  sin  oj  -  sin  Cl  cos  oj  cos  i  sin  Cl  sin  i 
sin  Cl  cos oj  +  cos  Cl  sin  oj  cos  i  -sin  Cl  sinw  4-  cos  Cl  cos  oj  cos  i  -cos  Cl  sin  i 
sin  oj  sin  i  cos  oj  sin  i  cos  i 


(B.2) 


For  notational  convenience  in  the  following,  vectors  with  no  superscript  will  be  assumed 
to  be  in  the  PQW  coordinate  system.  Those  vectors  which  are  written  with  only  two 
components  have  their  W-component  equal  to  zero;  that  is,  they  lie  in  the  orbital  plane. 

Two  derived  quantities,  used  in  computing  the  orbit,  are  the  mean  anomaly  at  time 
tT 


M  4  M0  +  vV*3  (t  -  t0)+27rk  ,  (B.3) 


t  From  a  developmental  viewpoint,  the  basic  quantity  is  u,  the  angle  subtended  at  the  focus  of  the 
ellipse.  The  mean  anomaly  M,  as  well  as  the  eccentric  anomaly  E,  are  derived  from  u. 


B-l 


where  0  <  M  <  2ir,  and  the  eccentric  anomaly  E  is  determined  implicitly  by  the  equa¬ 
tion 


E  -  e  sin  E  =  M  .  (B.4) 

The  constant  /<  is  the  gravitational  parameter  (gravitational  constant  times  the  mass  of 
the  earth,  (lj),  and  the  quantity  'J \if  a3  is  termed  the  mean  motion  since  it  is  the  rate  of 
change  of  the  mean  anomaly  M. 

The  position  of  the  satellite  at  time  t  is  given  by 

r  =  (r  cos  u,  r  sin  u)  ,  (B.5) 

where  r,  the  magnitude  of  r,  is  determined  from  ([1],  equation  (4-2-14)) 

r=  ||  r  ||  =  a  (1  -  e  cos  E) ,  (B.6a) 


and  v,  the  angle  subtended  by  r  with  the  P-axis  at  the  focus  of  the  ellipse,  satisfies  (|lj, 
equations  (4-2-12)  and  (4-2-13)) 


cos  v  —  —  (cos  E  -  e) 
r 

(B.6b) 

sin  v  =  —  \/l  -  e2  sin  E  . 
r 

(B.6c) 

Thus, 

r  =  a  (cos  E  -  e,  \/l  -  e2  sin  E) 

(B.7) 

Then, 

r  =  (-E  a  sin  E,  E  a  \/l  -  e2  cos  E)  . 

(B.8) 

But,  differentiating  (B.4)  with  respect  to  t,  using  M  ■=  Vfi/a?  from  (B.3), 
ship  (B.6a),  we  have 

and  relation- 

e=>/£, 

r  V  a 

(B.9) 

so  that 


B-2 


r  =  (-sin  E,  \/l-e2cosE) 
r 


(B.10) 


or 


-sin  E 


cos  E 


,  Vi¬ 


cos  E 


1  -  e  cosE 


(B.11) 


Next  we  compute  the  partial  derivatives  of  r  and  r  with  respect  to  the  elements. 
Equation  (B.4)  implies 


i.e., 


dE  (1  -  e  cos  E)  -  sin  E  de  =  dM  , 


1 


dE  = _ _  _ 

dM  1  -  e  cos  E  r 


(B.12) 


(B.13a) 


dE  = _ 

de  1  -  e  cos  E 


sin  E  a  sin  E 


(B.13b) 


For  convenience,  let 


m 


(B.14) 


be  the  elapsed  mean  motion.  Considering  M  a  function  of  a  and  t,  we  have 


dM 

SMn 


=  1, 


(B.15) 


and 


From  (B.13a) 


dM 

da 


-}V2F(t-g-lt 

dM  _  nr 
dt  V  a3  ’ 


dE  dE  dM  _  a  3  m  _  3  m 

da  dM  da  r  2  a  2  r 


(B.1G) 


(B.X7) 


(B.18) 


We  may  now  compute  the  partials  of  r  with  respect  to  the  elements  a,  e,  M0. 


B-3 


=  (cos  E  -  e,  Vl  -  e2  sin  E)  +  a  (-sin  E,  Vl  -  e2  cos  E)  ~~ 

O  ct  0  3/ 


=  (cos  E  -  e,  \/l  -  e2  sin  E)  -  pm  (  — S‘n  ^  ,  Vl  -  e2  — co-  —  — 
v  '  2  y  1  -  e  cosE  1  -  e  cos  E 


cos 


17  *  3  sin  E  /-  2  ■  T?i  3  a  2  cos  E  |  /r>  1 0\ 

E+—  m- - — ,  v  1  -  e*  sm  E  -  — mv  1  -  e  - —  .  (B.19) 

2  1-ecosE  2  1  -  e  cos  E  J  ’ 


dr  {  .  -e  sin  E 

3e  (  Vl- e2 


—  +  a  (-sin  E,  Vl  -  e2  cos  E)-^- 
O  '  ’  '  de 


I  i  •  r,  sin  E  e  sin  E  ,  .  2  ■  ™  cos  E 

—  a  1-1-smE - — , - * - !-  Vl  -e  sin  E - 


1  -  e  cos  E  ’ 


\/l  -  e2 


1  -  e  cos  E 


(B.20) 


dr 

3M0 


=  a  (-sin  E, 


Vl  -  e2  cos  E) 


dE 

dM0 


-sin  E  y/j  _  e2  cos  E 
1  -  e  cos  E  ’  1  -  e  cos  E 


(B.21) 


since  3M/3M0  =  1. 

Next,  we  compute  the  derivatives  of  the  velocity  r.  First,  from  (B.10), 


dr  _  r  d r 

dE  ~  r  5E 


-^(-cosE,  -Vl  -  e2 

r 


sin  E) 


or,  by  (B.6a), 

dr  _ e  sin  E  //T|  cos  E  p'j  ^2  sin  E 

dE  1  -  e  cos  E  V  a  |  1-ecosE’  1  -  e  cos  E 


(B.22) 


(B.23) 


Note  that  equations  (B.22)  and  (B.23)  were  computed  by  considering  f  as  a  function  of 
a,  e  and  E.  That  shall  be  taken  to  be  the  meaning  of  dr/dE,  i.e.,  that  quantity  is  given 
by  (B.22)  or  (B.23).  However,  we  wish  to  consider  a,  e,  M0  as  the  independent  variables. 
Thus,  below,  considering  E  as  a  function  of  a  and  Mo,  we  make  use  of  3E/3a  as  given  by 
(B.18).  We  must  be  careful  to  be  consistent,  however,  since  f  is  not  expressed  uniquely 
in  terms  of  a,  e,  E,  and  M0.  The  route  used  for  the  chain  rule  will  be  equation  (B.10) 
with  r  from  equation  (B.6a)  followed  '  •  substitution  of  E  as  determined  from  (B.4)  and 
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From  (B.10),  (B.6a),  and  (B.18) 


(B.3) 


dr 

3a 


T_ 

r 


(1  -  e  cos  E)  + 


3r_3E 
dE  3a 


r  3  m  3r 

2a  2  r  3E 


if £,  3  m _ 3r 

a  ^  2  2  1  -  e  cos  E  3E 


Also, 


di- 

3e 


(-a  cos  E)  + 


-cos  E  |  3r  3E 

^  _  e2  6  +  3E  3e 
* 


(B.24) 


Lastly, 


cos  E 

1  -  e  cos  E 


f 


e _ cos  E 

\J\-  £  1  -  e  cos  E 


+ 


sin  E 

1  -  e  cos  E 


dr 

3E' 


(B.25) 


dr  _  dr  3E  __  1 _ 3r 

3M0  3E  3M0  1  -  e  cos  E  3E 


(B.26) 


To  compute  the  derivatives  of  r  and  r  in  UK  coordinates  with  respect  to  an  arbi¬ 
trary  element  q,  we  use  (B.l): 


3r  _ 3r 

3q(IJK)  ~~  5q(PQW) 


for  q  =  a  or  e  or  M0 


(B.27) 


3r  _  3U_  (pqw) 
3q^IJK)  3q 


for  q  =  i  or  fl 


(B.28) 


* 


d\J 

he  same  equations  hold  for  r.  It  remains  to  compute 

3q 


We  have 


3U 

3i 


sin  Cl  sin  w  sin  i  sin  Cl  cos  w  sin  i 
-cos  fl  sin  u>  sin  i  -cos  Cl  cos  u>  sin  i 
sin  o>  cos  i  cos  u>  cos  i 


sin  Cl  cos  i 
-cos  w  cos  i 
-sin  i 


(B.29) 


3U 

3f2 


-U2i  -U22  -U23 

Un  U12  u13 

0  0  0 
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(B.30) 


au 

du> 


U 

U 

U 


12 

-U,1 

22 

-Uo! 

32 

-U31 

0 

0 

-sin  i 


(B.31) 


B-6 


